AD-A040  843 


UNCLASSIFIED 


UNITED  TECHNOLOGIES  RESEARCH  CENTER  EAST  HARTFORD  CONN  F/G  20/4 

A STUDY  OF  THE  TURBULENT  SHOCK  WAVE  BOUNDARY  LAYER  INTERACTION. CU) 
FEB  77  R LEVY*  S J SHAMROTH.  H J GIBELING  F33615-73-C-4119 


UTRC/R76-91 1 729-15 


AFFDL-TR-76-163 


o 

Tf* 

© A STUDY  OF  THE  TURBULENT  SHOCK  WAVE  BOUNDARY 
< LAYER  INTERACTION 

CD 

«=c 


UNITED  TECHNOLOGIES  RESEARCH  CENTER 


TECHNICAL  REPORT  AFFDL  TR  76-163 

FINAL  REPORT  FOR  PERIOD  JULY  1,  1973  - OCTOBER  25,  1976 


0 

CJ> 

1 , 1 AIR  FORCE  FLIGHT  DYNAMICS  LABORATORY 


. -_J  AIR  FORCE  WRIGHT  AERONAUTICAL  LABORATORIES 

O L_i_  AIR  FORCE  SYSTEMS  COMMAND 
2S  WRIGHT  PATTERSON  AIR  FORCE  BASE,  OHIO  45433 


\ 

' NOTICE 

- 

When  Government  drawings,  specifications,  or  other  data  are  used  for 
any  purpose  other  them  in  connection  with  a definitely  related  Government 
procurement  operation,  the  United  States  Government  thereby  incurs  no 
responsibility  nor  any  obligation  whatsoever;  and  the  fact  that  the 
government  may  have  formulated,  furnished,  or  in  any  way  supplied  the 
said  drawings,  specifications,  or  other  data,  is  not  to  be  regarded  by 
implication  or  otherwise  as  in  any  manner  licensing  the  holder  or  any 
other  person  or  corporation,  or  conveying  any  rights  or  permission  to 
manufacture,  use,  or  sell  any  patented  invention  that  may  in  any  way  be 
H O related  thereto. 

i 

The  work  described  in  this  Final  Technical  Report  was  supported  by 
the  U.S.  Air  Force  through  Dr.  Wilbur  Hankey  of  the  Air  Force  Systems 
Command,  Wright- Ratterson  AFB,  Ohio,  under  Contract  F336l5-73-C-4ll9, 
i Project  70640307,  and  Work  Unit  23070419.  Dr.  Joseph  Shang  (AFFDL/FXM) 

i was  the  contract  monitor  and  the  work  covered  by  this  report  was  performed 

at  United  Technologies  Research  Center,  East  Hartford,  Connecticut, 
between  July  1,  1973  and  October  25,  1976.  The  final  report  was  submitted 
February  18,  1977- 


This  report  has  been  reviewed  by  the  Information  Office  (10)  and  is 
releasable  to  the  National  Technical  Information  Service  (NTIS).  At  NTIS, 
it  will  be  available  to  the  general  public,  including  foreign  nations. 

This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


A|k  »0RO 


< MAY  77  - 150 


RFAD  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


REPORT  DOCUMENTATION  PAGE 


t report  wciW'pm- 

AFFDL-7E-76-163 


■i  PECiPrn  JJi  CATA_OG  N UM  BE  R 


GOVT  ACCX.SA1CN  HO. 


- TVFT  OF  REPORT  & PERIOD  COVERED 

Final  Report  , 

July  1,  1973-October  25,  19' 


• TITLE  rand  SuOI/ll.;  " *" 

A Study  of  the  Turbulent  Shock.  Wave  Boundary  J 
Layer  Interaction, 


^~-P-S«fcOEU<J!tS.-0H0  REPORT  NUMBER 

R76-9H729-15 

T CONlfiktl  UPG R~A NT  HUMBERTS' 


7.  AUTHOR(a) 


Ralph/Levy,  Stephen  J./shamroth, 
Howard  J.l  Gibeling,  Henry^/tocDonald 


9 PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 


APE  A A WORK  UNIT  NUMBERS 

Project  706L0307  If 
Work  Unit  23070419 ; 


United  Technologies  Research  Center 
East  Hartford,  Connecticut  06ll8 


Air  Force  Flight  Dynamics  Laboratory  -/  V 

Air  Force  Systems  Command 

Wright-Patterson  Air  Force  Base,  Ohio  4^433 

• 4 MONITORING  AGENCY  NAME  h ADDRESS^//  different  from  Controlling  Office) 


OECL  ASSlFlC  ATION  DOWNGRADING 
SCHEDULE 


16  DISTRIBUTION  STATEMENT  (of  this  Report) 


Approved  for  Public  Release;  Distribution  Unlimited 


17.  DISTRIBUTION  STATEMENT  (of  (he  abstract  entered  in  Block  20.  If  different  from  Report) 


18  SUPPLEMENTARY  NOTES 


19  KEY  WOROS  ( Continue  on  reverse  side  if  necessary  and  identify  by  block  number) 


Viscous- Inviscid  Interaction,  Shock  Wave-Boundary  Layer  Interaction 
Navier-Stokes  Solutions . 


20  ABSTRACT  ( Continue  on  reverse  side  It  necessary  and  identify  by  block  number) 


Two  approaches  to  the  shock  wave-boundary  layer  interaction  procedure 
have  been  investigated.  In  the  first  approach  a strong  interaction  boundary 
layer  procedure  has  been  applied  to  the  shock  wave -turbulent  boundary  layer 
interaction  problem  by  modifying  a well-proven  weak  interaction  boundary  layer 
procedure  to  allow  for  strong  inueraction  solutions.  -sThe  modifications 
included  the  incorporation  of  an  outer  flow  interaction  model  and  the  revision 


FORM 

I JAN  1! 


UNCLASSIFIED 


EDITION  OF  I NOV  65  IS  OBSOLETE 


ECUHITl  CLASSIFICATION  OF  THIS  P AGE(IFh*n  O.l.  Enffd) 


of  the  governing  equations  to  allow  for  flow  in  separated  regions.  When  a 
turbulence  model  originally  developed  for  attached  boundary  layers  was 
used,  poor  agreement  was  found  between  calculations  and  data.  However, 
when  a modified  turbulence  model,  based  upon  computational  turbulence 
studies  and  recent  experimental  data  was  used,  considerably  better  agree- 
ment was  obtained.  The  original  weak  interaction  solution  procedure, 
the  modifications  made  under  the  current  effort,  and  sample  calculations 
are  presented  in  this  report.  'In  the  second  approach  the  feasibility  of 
using  a fully  implicit  Navier-Stokes  solution  procedure  for  the  inter- 
action problem  was  investigated.  The  results  indicate  that  the  fully 
implicit  procedure  is  a promising  one,  however,  further  development  work, 
particularly  in  regard  to  shock  representation,  is  required. 


^FCUWlTY  CLASSIFICATION  OF  this  PAOF'IIT'*'’  Fnf*r#<fl 


TABLE  OF  CONTENTS 


Section  Page 

T.  INTRODUCTION 1 

II.  DISCUSSION 6 

THE  BOUNDARY  I AYER  APPROACH 6 

Weak  Interaction  Boundary  Layer  Procedure 6 

Modifications  For  Strong  Interaction  9 

Solution  Of  The  Strong  Interaction  Boundary  Layer 

Equations i4 

Results 17 

THE  NAVIER-STOKES  APPROACH 41 

The  Basic  Analysis 42 

Method  of  Solution 47 

Results  of  Navier -Stokes  Computations 48 

III.  CONCLUSIONS 54 

BOUNDARY  LAYER  ANALYSIS  B4 

NAVIER-STOKES  ANALYSIS 54 

APPENDIX  A:  WEAK  INTERACTION  CALCULATION  PROCEDURE  55 

APPENDIX  B:  BOUNDARY  LAYER  MODIFICATIONS  FOR  6*  SPECIFIED 57 

APPENDIX  C:  ORIGINAL  TURBULENCE  MODEL  USING  THE  TURBULENCE 

KINETIC  ENERGY  EQUATION  5° 

APPENDIX  D:  THE  MINT  CALCULATION  PROCEDURE 63 


J 


I.IST  OF  ILLUSTRATIONS 


Page 

Departure  Solutions  to  the  Boundary  Layer  Equations  ....  15 

Experimental  Free  Interaction  Pressure  Rise  19 

Laminar  Free  Interaction  Calculations  21 

Laminar  Calculation  Branch  Control  Using  Kick  22 

Laminar  Calculation  Branch  Control  Using  X^oRNER  23 

Laminar  Strong  Interaction  Flow  Over  Compression  Corner  . . 24 

Separating  Turbulent  Strong  Interaction  Solution  Using 

Weak  Interaction  Turbulence  Modeling  Assumptions  26 

Comparison  of  Equilibrium  and  TKE  Options  For  Attached 

Boundary  Layer  Turbulence  Model  27 

Wake  Mixing  Length  Distribution  with  Attached 
Boundary  Layer  Turbulence  Model  Including  Low 

Reynolds  Number  Correction  29 

Sensitivity  of  Free  Interaction  to  Turbulence 

Modeling  Assumptions  30 

Mixing  Length  Distributions  with  Attached  Boundary 

Layer  Turbulence  Model 31 

Free  Interaction  Wall  Pressure  Distributions  33 

"Equilibrium"  Wake  Mixing  Length  Distributions 

with  Low  Reynolds  Number  Correction 34 

"Equilibrium"  Normal  Mixing  Length  Distributions 

with  Low  Reynolds  Number  Corrections  35 

"TKE"  Wake  Mixing  Length  Distributions  with  Low 

Reynolds  Number  Correction  36 

"TKE"  Normal  Mixing  Length  Distributions  with  Low 

Reynolds  Number  Corrections  37 

Turbulent  Strong  Interaction  Solution  For  19.67°  Ramp 

Using  Modified  Turbulence  Model  36 

Turbulent  Strong  Interaction  Solution  For  l6.06°  Ramp 

Using  Modified  Turbulence  Model  39 

Turbulent  Strong  Interaction  Solution  For  9*81°  Ramp 

Using  Modified  Turbulence  Model  4o 

Navier-Stokes  Interaction  Calculation,  Pg/P^  = 2 52 

Navier-Stokes  Interaction  Calculation,  ^2/^1  =3  53 


SECTION  I 
INTRODUCTION 

Viscous-inviscid  interaction  flow  fields  are  common  occurrences  found 
both  on  high  Mach  number  flight  vehicles  and  in  supersonic  or  transonic  turbo- 
machinery. in  general,  the  interaction  may  result  either  from  the  impingement 
of  a shock  wave  upon  a boundary  layer  or  from  a boundary  layer  developing 
along  a compression  surface.  The  resulting  flow  field  is  a complex  phenomenon 
governed  by  the  mutual  interaction  of  a viscous  inner  shear  layer  and  a 
nominally  inviscid  outer  flow.  Since  viscous-inviscid  interactions  may  lead 
to  high  drag  losses,  loss  of  control  effectiveness,  and  very  high  local 
heating  rates,  and  since  such  interactions  commonly  occur  on  high  Mach  number 
flight  vehicles,  a method  to  predict  the  interaction  flow  field  accurately 
would  be  a valuable  asset  to  the  vehicle  design  team.  The  interaction  flow 
field  occurring  in  most  cases  of  practical  interest  involves  a turbulent 
boundary  layer  and,  thus,  a calculation  procedure  capable  of  accurately 
modeling  turbulent  boundary  layer  development  should  form  the  basis  of  the 
interaction  prediction  procedure. 

Following  Prandtl, the  now  classical  method  of  solving  viscous  boundary 
layer-type  problems  first  calculates  a streamwise  pressure  distribution  from 
purely  inviscid  considerations  and  then  calculates  the  boundary  layer  develop- 
ment under  the  influence  of  this  inviscid  pressure  distribution.  No  attempt 

is  made  *o  correct  t lie  inviscid  pressure  calculation  for  viscous  effects, 
originating  in  the  boundary  layer  region.  Obviously,  such  a two-step  pro- 
cedure can  be  valid  only  if  the  pressure  distribution  calculated  ignoring 
viscous  effects  is  a sufficiently  good  approximation  to  the  pressure  dis- 
tribution which  actually  occurs.  Tn  many  cases  of  practical  interest,  ‘he 
inviscid  pressure  distribution  is  an  excellent  approximation  to  physical 
reality  and  in  these  cases  the  classical  two-step  procedure  is  valid.  How- 
ever, in  certain  flow  situations  which  arc  termed  strong  interaction  flows, 
the  viscous  displacement  effects  near  the  body  surface  are  sufficiently  large 
'o  cause  the  pressure  distribution  calculated  ignoring  viscous  wall  offer's 
to  he  in  serious  error  and  in  these  cases  the  classical  two-step  boundary 
layer  calculation  procedure  is  invalid.  The  shock  wave  boxindary  layer 
interaction  represents  a common  type  of  strong  interaction  flow.  For  such 
a flow,  a prediction  procedure  which  can  be  used  wi‘h  confidence  must  accoun* 
for  'he  mutual  in' enaction  between  the  inner  viscous  flow  and  the  outer 
s-innlly  inviscid  flow. 

Most  interaction  flow  field  procedures,  which  have  been  developed  and  are 
discussed  in  the  open  literature  have  extended  classical  boundary  layer 
theory  to  account  for  the  mutual  interaction  phenomena.  Under  this 
approach  the  boundary  layer  equations  which  represent  * he  inner  flow  fie] i 
are  solved  simultaneously  with  an  equation  representing  the  outer  flow  field 
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thus  including  the  required  mutual  interaction  in  the  solution  process. 

Most  strong  interaction  solutions  of  the  boundary  layer  equations  are  based 
conceptually  upon  the  original  work  of  Crocco  and  Lees  (Ref.  13)  as  modi- 
fied by  Lees  and  Reeves  (Ref.  °7 ) . In  brief,  Lees  and  Reeves  solved  a set 
of  integral  laminar  boundary  layer  equations  under  the  constraint,  that  the 
flow  angle  at  the  outer  edge  of  the  boundary  layer  as  determined  by  the 
boundary  layer  equations  be  equal  to  the  flow  angle  at  the  outer  edge  of  the 
boundary  layer  as  determined  by  the  inviscid  flow  field  equations.  This 
constraint  was  expressed  through  an  equation  relating  the  outer  edge  flow 
angle  as  determined  by  the  continuity  equation  to  ^he  angle  determined  by 
the  Prandtl-Meyer  law.  The  use  of  an  additional  equation  allowed  the  stream- 
wise  static  pressure  distribution  to  be  calculated  simultaneously  with  the 
boundary  layer  development  and,  thus,  the  streamwise  static  pressure  dis- 
tribution emerged  as  part  of  the  solution  as  it  should.  The  Lees  and  Reeves 
procedure  resulted  in  fairly  good  predictions  of  static  pressure  through  the 
int eract ion  region  for  laminar  boundary  layer  development  on  adiabatic  walls; 
however,  Use  prediction  of  skin  friction,  although  qualitatively  correct, 
showed  significant  quantitative  disagreement  with  experimental  data.  The 
Lees  and  Reeves  integral  theory  was  ext. ended  to  nonadiabatic  walls  by 
Klineberg  and  Lees  (Ref.  24)  and  Holden  (Ref.  22). 

Although  the  analyses  of  Refs.  27,  22  and  24  have  proven  capable  of 
predicting  the  streamwise  pressure  distribution  through  the  interaction 
region,  they  are  all  based  upon  integral  solutions  of  the  equations  and, 

' hus,  are  limited  in  their  generality.  A more  general  solution  for  the 
interaction  problem  has  been  developed  b>  Reyhner  and  Flugge-Lotz  (Ref.  4l) 
who  replaced  the  integral  boundary  layer  equations  with  the  boundary  layer 
partial  differential  equations  of  motion  thus  eliminating  the  constraints 
imposed  by  the  profile  assumptions  required  by  integral  solutions.  Since 
the  boundary  layer  equations  calculated  by  spatial  marching  procedures  are 
unstable  in  regions  of  reversed  flow,  the  finite  difference  algorithm  was 
modified  so  that  in  the  separated  flow  region  the  streamwise  convective 
terms  were  ignored  or  given  a false  positive  convection  velocity  to  main- 
tain numerical  stability.  Being  a full  finite-difference  solution  rather 
than  an  integral  solution,  the  Reyhner  and  Flugge-Lotz  procedure  is  more 
general  than  integral  calculation  procedures;  however,  like  integral  pro- 
cedures, the  Reyhner-Flugge-Lotz  procedure  is  still  confined  to  laminar 
flow. 


A solution  of  the  turbulent  boundary  layer  equations  for  strong  inter- 
actions has  been  presented  by  Bertke,  Werle  and  Polak  (Ref.  5 ) who  solved  a 
set  of  equations  originally  developed  by  Werle  and  Vatsa  (Ref.  57 ) for  laminar 
flows.  The  method  of  Ref.  77  is  apseudo-time  relaxation  scheme  that  allows  adirect 
specification  of  the  downstream  boundary  conditions  without  having  to  adjust 
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initial  parameters.  A turbulence  model  that  had  been  developed  for  attached 
weak  interaction  boundary  layer  flows  was  used.  When  this  method  was  applied 
to  the  turbulent  free  interaction  problem,  the  calculated  plateau  pressures 
were  found  to  be  approximately  30  percent  higher  than  available  data  correla- 
tions. Bertke,  Werle,  and  Polak  attributed  the  discrepancy  in  predicted 
plateau  pressure  to  limitations  of  the  turbulence  model  used.  A broad  review 
of  shock  wave-boundary  layer  interaction  solutions  has  been  compiled  byHankey  (Ref.  21 ). 

Although  until  recently  numerical  viseous-inviscid  interaction 
calculations  were  based  almost  entirely  upon  ex4 ended  boundary  layer  pro- 
cedures, 4 he  rapid  development  of  both  numerical  techniques  and  corrrpu4  er 
capability  now  have  made  it  practical  to  apply  Navier-Stokes  calculation 
procedures  to  the  interaction  problem.  Since  4he  Navier-Stokes  equations 
represent  the  exact  equations  of  mo* ion,  they  contain  no  approximations 
o4her  than  those  made  in  connection  with  *urbulence  modeling.  For  example, 
the  assumption  of  constant  pressure  at  any  streamwise  sta*ion  which  is 
inheren4  in  most  boundary  layer  procedures  and  which  may  be  in  serious 
error  in  supersonic  flow  is  relieved  in  the  Navier-St.okes  approach.  Further- 
more, the  approximat ion  of  the  convective  terms  in  separated  flow  regions 
which  is  necessary  for  numerical  stability  in  boundary  layer  calculat ions 
is  no4  required  in  Navier-Stokes  calculations.  Thus,  the  basic  equations 
used  *o  solve  the  interaction  problem,  via  a Navier-Stokes  calcula4 ion  are 
not  compromised  by  approximations  which  are  required  for  a successful 
boundary  layer  solution. 

The  advan* ages  of  the  Navier-Stokes  solution  must  be  balanced  against 
cer4 ain  disadvantages  such  as  code  complexity,  code  run  4 ir.e  and  code 
storage  requiremen* s.  Navier-Stokes  calculation  procedures  in  general 
require  more  code  development  time  than  do  boundary  layer  procedures.  In 
addition,  Navier-Stokes  procedures  may  require  increased  computer  storage 
and  increased  computer  run  time  to  ob+ain  a successful  solution.  A consid- 
eration of  both  the  advantages  and  disadvantages  of  Navier-Stokes  solutions 
indicate  that  the  potentially  increased  accuracy  of  this  approach  may 
well  balance  out  the  increased  requirements  in  computer  storage  and  run 
time.  With  these  considerations  in  mind  i*  seems  reasonable  to  continue 
development  of  both  approaches  to  the  interaction  problem  at  the  present 
time . 

Solutions  of  the  Navier-Stokes  equations  for  shock  wave- turbulent 
boundary  layer  flow  interactions  have  been  presented  recen* ly  by  a variety 
of  authors  including  Wilcox  (Ref.  56),  Baldwin  and  MacCormack  (Ref.  ?),  and 
Shang  and  Hankey  (Ref.  !(7).  In  Ref.  56,  Wilcox  combined  an  explicit 
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first-order  finite-difference  scheme  with  Saffman's  turbulence  model 
(Ref.  UU ) to  predict  two-dimensional  turbulent,  supersonic  flows  with 
separation.  Although  Wilcox  was  successful  in  obtaining  flow  field  solu- 
tions, the  predicted  surface  pressure  distribution  disagreed  significantly 
with  experimen* al  data.  In  a more  recent  work,  Baldwin  and  MacCormack 
(Ref.  .?)  applied  MacCormack's  two-step  second-order  explicit  Navier-St OKe. 
method  4o  the  shock  wave-turbulent,  boundary  layer  interaction  problem. 
Although  both  a mixing  lerpgt  h turbulence  model  and  a 'wo-equation  trans- 
port model  of  turbulence  were  used,  neither  turbulence  model  produced  flow 
fields  that  agreed  well  with  measurements.  Shang  and  Hankey  applied  4 he 
explicit  two-dimensional  time-dependent  method  of  MacCormack  to  the  inter- 
act ion  flow  field  problem  using  an  eddy  viscosity  model  which  allowed  for  a 
la:  in  the  response  of  Hie  4urbulence  to  the  sudden,  severe  adverse  pressure 
gradient  tha4  is  characteristic  of  strong  interacting  separating  boundary 
layers.  A rate  relation  was  introduced  to  decay  the  eddy  viscosity  expon- 
entially from  the  value  at  the  start  of  the  interaction  toward  the  local 
equilibrium  value  with  a time  scale  correlated  *o  the  boundary  layer  t hick- 
ness. Utilising  the  relaxation  model  for  eddy  viscosity,  Chan:  and  Hankey 
obtained  predictions  of  a series  of  compression  corner  flows  which  were  in 
good  agreement  with  the  experimental  data  of  Law  (Ref.  16).  More  recently, 
Chang,  Hankey  and  law  (Ref.  H8)  have  compared  predictions  for  a series  of 
ncident  shock  wave-boundary  layer  interactions  wi'  h experimental  da4 a. 
Again,  good  agreement  was  obtained  between  the  calculations  and  4 he  experi- 
ment . 


The  Navier-Stokes  solut ions  discussed  in  t he  previous  paragraph  were  all 
based  upon  an  explicit  difference  scheme  for  solving  the  unsteady  form  of 
the  governing  equations.  An  initial  flow  field  is  assumed  and  the  calcula- 
tion ‘hen  proceeds  in  time  until  a steady  state  is  reached.  However,  explicit 
numerical  procedures  of  this  type  are  subject  to  one  or  more  stability 
restrictions  on  the  size  of  the  time  step  relative  to  the  spatial  mesh  sine. 
These  stability  limits  usually  correspond  to  the  well  known  Courant- 
Fredricks-Lewy  (CFL)  condition  and  in  some  schemes  to  an  additional  stability 
condition  arising  from  viscous  terms.  These  stability  restrictions  can  lower 
computational  efficiency  by  imposing  a smaller  time  step  than  would  be 
otherwise  desirable.  Since  the  restriction  becomes  more  and  more  stringent 
as  the  spatial  :rid  is  refined,  the  restriction  can  become  particularly 
burdensome  in  the  calculation  of  a turbulent  boundary  layer  where  a very 
fine  mesh  near  the  wall  may  be  required.  In  contrast  to  most  explicit- 
methods,  implicit  methods  tend  *o  be  stable  for  large  4 ime  steps  and 
therefore  offer  the  prospect  of  substantial  increases  in  comput at ional 
efficiency.  When  the  present  work  was  initiated,  no  implicit  Navier-Stokes 
procedure  for  nredict  itig  1 lie  shock  wave-boundary  layer  interact  ion  flow 
field  had  been  reported  in  the  open  literature  and  in  view  of  the  potential 
benefits  associated  with  implicit  methods,  il  was  decided  *o  apply  an 


4 


1 


I 

e 


i 

I 


implicit 
a t tidy  1 o 
realized 
could  bo 
implici* 


technique  to  this  problem.  This  application  represents  a feasibility 
determine  if  the  potential  benefits  of  implicit  methods  could  be 
in  an  interaction  problem.  If  it  appears  that  substantial  benefits 
realized,  then  future  work  could  tie  aimed  4 oward  developin'-'  the 
interaction  solution  procedure  as  a design  tool. 


Since  the  initiation  of  the  present  work, a technique  in  which  viscous 
terms  are  treated  implicitly  and  convective  terms  are  treated  explicitly 
has  been  presented  by  MacCormack  (Ref.  30).  The  present  Navier-Stokes 
investigation  centers  on  a fully  implicit  procedure,  the  Multi- dimensional 
Implicit.  Nonlinear  Time-dependent  (MINT)  technique,  developed  by  Briley 
and  McDonald  (Ref.  b) . The  method  has  been  ex4  ended  to  turbulent  flow  by 
Briley,  McDonald  and  Gibeling  (Ref.  10 ) and  to  flow  in  combustors  by 
Gibeling,  McDonald  and  Briley  (Ref.  18).  In  the  Navier-Stokes  portion 
of  the  present  effort  a two-dimensional  version  of  the  MINT  procedure 
is  applied  4o  the  shock  wave-boundary  layer  interaction  problem. 


The  present  report  covers  all  work  performed  under  the  subject  contrac4- 
and  can  be  divided  into  two  parts.  The  first  portion  concerns  the  boundary 
layer  approach;  this  represents  the  major  portion  of  the  effort  done  under 
the  subject  contract.  Under  the  boundary  layer  approach,  a well-proven  weak 
interaction  boundary  layer  procedure  was  modified  to  allow  for  strong  inter- 
action solutions.  The  modifications  include  (i)  the  incorporation  of  an  outer 
flow  analysis,  (ii)  the  revision  of  the  governing  equations  to  allow  for  flow 
in  separated  regions,  and  (iii)  the  revision  of  the  turbulence  model.  The 
original  weak  interaction  solution  procedure,  the  modifications  made  under 
the  current  effort,  and  sample  calculations  are  described  in  the  following 
sections.  The  second  portion  of  the  report  concerns  the  Navier-Stokes 
approach.  The  governing  equations  and  the  numerical  method  are  discussed  and 
sample  calculations  are  presented. 
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SECTION  II 
DISCUSSION 


ME  BOUNDARY  LAYER  APPROACH 

’’he  asis  for  the  strong  i iteractiou  boundary  layer  calculations  described 
in  this  repor * is  a well-proven  weak  interaction  boundary  layer  procedure  which 
has  been  used  to  compute  successfully  a wide  variety  of  laminar,  transitional, 
and  turbulent  boundary  layers.  A complete,  detailed  description  of  the  procedure 
as  well  as  an  extensive  comparison  between  predictions  and  iata  can  be  found  in 
fs.  25,  33*  and  46.  The  present  discussion  first  briefly  describes  the 
weak  interaction  procedure  and  then  proceeds  to  detail  the  modifications  re- 
quired for  a strong  interaction  code.  These  modifications  include  the  incor- 
pora' ion  of  an  outer  flow  analysis  into  the  calculation,  the  modification  of 
'he  equations  in  the  reversed  flow  region,  and  the  modification  of  the  ‘urbu- 
lence  model  within  the  separation  bubble.  Next,  a discussion  of  the  cause  and 
characteristics  of  the  so-called  branching  solutions  obtained  with  the  strong 
interaction  code  is  presented.  The  validity  of  the  branching  solutions  along 
with  the  characteristics  of  free  interactions  are  presented  an!  the  behavior  of  the 
branches  after  encountering  a corner  or  an  incident  shock  wave  are  demonstrated. 
Next,  attention  is  focused  upon  the  downstream  boundary  condi: ion  problem.  The 
boundary  layer  calculation  approach  used  to  solve  the  strong  interaction  problem 
requires  he  downstream  boundary  conditions  to  be  satisfied  through  an  iterative 
'shooting'  procedure.  In  this  iterative  procedure  a sequence  of  initial  value 
problems  is  solved  until  the  solution  at  the  end  of  the  strong  interaction  regioii 
is  compatible  with  specified  downstream  boundary  conditions.  The  initial  value 
problems  within  the  sequence  differ  in  the  value  of  a free  parameter.  Different 
values  of  the  free  parameter  lead  to  different  branching  solutions  and  ulti- 
mately to  different  flow  fields  at  the  downstream  flow  field  boundary.  The 
choice  of  the  free  parameter  is  discussed  and  the  rationale  for  using  one 
variable  for  the  parameter  when  the  solution  is  far  from  convergence  and 
another  variable  for  the  parameter  when  the  solution  is  near  convergence  is 
presented.  Finally,  the  results  of  applying  the  methodology  outlined  above 
are  presented.  Both  laminar  and  turbulent  results  are  presented  along  with 
the  details  of  the  modified  turbulence  model  that  was  used  in  the  successful 
turbulent  strong  interaction  calculations. 

Weak  Interaction  Boundary  Layer  Procedure 

The  weak  interaction  solution,  upon  which  the  current  work  is  based,  has 
been  developed  over  a period  of  years  and  details  of  the  procedure  are  presented 
tp  Refs.  25,  33,  34  and  46.  The  procedure  itself  solves  the  steady  state 
boundary  layer  equations  for  laminar,  transitional,  and  turbulent  flow.  As 
shown  by  -any  authors  (for  example,  Schubauer  and  Tchen  (Ref.  45)),  for  two- 
dimensional  or  axisymnei ric  flows,  steady  in  the  mean,  the  bovmdary  layer 
approximations  to  the  momentum,  energy,  and  continuity  equation  become: 
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In  Eqs.  (l)  through  (3)  x and  y are  coordinates  in  the  streamwise  and  trans- 
verse directions,  u and  v are  velocity  components  in  the  x and  y directions, 

0 is  density,  P is  pressure,  C is  specific  heat,  T°  is  total  temperature,  r 
is  the  radius  of  curvature  for*  an  axisymmetric  body,  and  the  exponent  a is  zero 
for  two-dimensional  flows  and  unity  for  axisymmetric  flows.  By  definition,  the 
effective  shear  stress,  t,  and  the  effective  heat  transfer,  Q are  given  as: 
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where  v is  viscosity,  k is  thermal  conductivity,  T is  static  temperature. 

In  Eqs.  (l)  through  (5)  overbars  indicate  averaged  quantities  and  primes  in- 
dicate fluctuating  quantities.  The  equations  are  valid  for  laminar,  transi- 
tional, or  turbulent  flows;  obviously,  for  laminar  flows  the  primed  quantities 
are  zero.  In  the  case  of  turbulent  and  transitional  boundary  layers,  it  is 
convenient  to  represent  the  contribution  of  the  apparent  turbulent  stress,  t^, 
to  the  total  shear  stress,  t,  by  an  effective  turbulent  viscosity,  i/j.  In  an 
analogous  manner,  the  turbulent  contribution  to  the  total  heat  flux,  Q is 
represented  by  an  effective  turbulent  conductivity,  krp,  such  that 
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and  by  analogy  to  laminar  flow  a turbulent  Prandtl  number,  PrT?  can  be  defined 

as 


At  this  point,  i/^  and  k,„  simply  represent,  definitions  and  in  no  way  limit  the 
subsequent  turbulence  modeling.  When  Eqs.  (h)  through  (7)  arc  substituted  into 
Eqs.  (l)  and  (?),  the  resulting  boundary  layer  momentum  and  energy  equal  ions 
take  the  form: 
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In  deriving  Lq.  (lO)  use  has  “been  made  of  the  definition  of  total  temperature 
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In  situations  in  which  the  flow  is  laminar,  Eqs.  (3)?  (9)>  and  (lO)  are 
solved  with  v,p  = 0 to  determine  the  boundary  layer  development.  When  the  flow 
is  transitional  or  turbulent,  it  is  necessary  to  model  vT  and  Pr^.  The  speci- 
fication of  the  turbulent  viscosity,  vT,  and  the  turbulent  Prandtl  number,  Prp 
is  carried  out  through  the  turbulence  model  described  in  detail  subsequently. 


In  the  present  procedure  the  equations  are  solved  by  using  the  continuity 
equation  to  eliminate  the  explicit  appearance  of  pv  from  the  momentum  and  energy 
equations.  Then  in  the  weak  interaction  case  where  the  streamwise  pressure 
gradient,  P(x)  is  specified,  the  momentum  and  energy  equations,  in  conjunction 
with  an  equation  of  state  and  equations  governing  and  Pr^,  form  a closed  set 
of  nonlinear,  parabolic,  partial  differential  equations  which  can  be  solved 
upon  specification  of  boundary  conditions.  Obviously,  in  the  strong  inter- 
action case  where  the  streamwise  pressure  gradient,  P(x)  is  not  specified 
a priori  but  must  emerge  from  the  solution,  an  additional  relation  is  required 
to  determine  the  static  pressure  distribution. 


The  wall  and  free-stream  boundary  conditions  employed  in  the  solution  are 
given  by: 
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The  subscripts  'w'  and  'e'  denote  wall  and  free-stream  conditions,  respectively. 
The  initial  conditions  for  the  problem  are  set  by  specifying  the  initial  boundary 
layer  displacement  thickness,  6*,  and  assuming  the  initial  development  is  sim- 
ilar in  the  dimensionless  coordinate,  T|,  where 

: ~ (lU) 

The  numerical  procedure  used  to  solve  the  ,'erning  momentum  and  energy 
equations  is  a Hartree-Womersley  approach  in  which  streamwise  derivatives  are 
replaced  by  finite  differences,  the  coordinate  normal  to  the  wall  is  nondim- 
ensionalized,  and  a stream  function  introduced.  The  resulting  momentum  equation 
is  a third  order  nonlinear  ordinary  differential  equation  in  the  transverse 
coordinate  and  the  energy  equation  is  a second  order  nonlinear  ordinary  differ- 
ential equation  in  the  same  transverse  coordinate.  At  each  streamwise  station 
the  nonlinear  coefficients  of  each  equation  are  estimated  from  the  solution  at 
the  previous  station  and  the  resulting  linearized  equations  solved  as  two-poinr- 
boundary  value  problems.  Appendix  A describes  the  solution  procedure  in  detail. 

The  fully- developed  turbulence  model  used  in  the  weak  interaction  procedure 
was  originally  presented  by  McDonald  and  Camarata  (Ref.  32)  for  two-dimensional 
incompressible  flow  and  has  been  extended  to  a variety  of  flow  situations  in 
Refs.  25,  33,  34  and  46.  The  model  is  based  upon  a solution  of  the  turbulence 
kinetic  energy  equation  which  is  a conservation  equation  derived  from  the  Navier- 
Stokes  equations  by  writing  the  instantaneous  quantities  as  a sum  of  mean  and 
fluctuating  parts.  The  i^h  Navier-Stokes  momentum  conservation  equation 
(i=  1,2,3,  referring  to  the  three  coordinate  directions)  is  multiplied  by  the 
i^h  component  of  fluctuation  velocity  and  the  average  of  the  resulting  three 
equations  is  taken.  The  three  averaged  equations  are  summed  to  obtain  the 
turbulence  kinetic  energy  equation.  The  derivation  of  the  turbulence  kinetic 
energy  equation  has  been  given  by  Favre  (Ref.  17)  for  compressible  fiow  and 
approximated  by  Bradshaw  and  Ferris  (Ref.  8)  to  boundary  layer  flows.  Two 
turbulence  models  are  used  in  the  current  strong  interaction  calculations;  the 
first  model  is  a modification  of  a Prandtl  mixing  length  equilibrium  turbulence 
model.  The  second  model  is  a modification  of  the  turbulence  kinetic  energy 
model  presented  in  Refs.  25,  33,  34  and  46.  The  modifications  are  discussed 
subsequently. 

Modifications  For  Strong  Interaction 

Matching  and  Marching  Procedures 

Although  the  weak  interaction  boundary  layer  procedure  described  above 
has  successfully  predicted  the  development  of  a wide  variety  of  boundary  layers, 
it  obviously  requires  as  input  a specified  pressure  distribution,  P(x).  In 
weak  interaction  flows  the  pressure  gradient  can  be  obtained  from  an  inviscid 
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calculation  which  ignores  all  viscous  displacement  effects.  However,  in  strong 
interaction  problems,  such  as  the  shock  wave  boundary  layer  interaction,  the 
inviscid  pressure  distribution  will  be  in  significant  error  as  the  pressure 
distribution  and  boundary  layer  development  mutually  affect  each  other.  The 
prediction  of  the  boundary  layer  development  and  the  pressure  distribution  using 
a procedure  which  recognizes  their  mutual  effects  upon  each  other  is  the  core 
of  the  basic  strong  interaction  problem. 

The  present  effort  follows  the  strong  interaction  approaches  of  Refs.  14, 
21,  24,  27  and  4l  which  introduce  the  additional  condition  tha  . the  boundary 
layer  flow  angle  must  agree  with  the  inviscid  flow  angle  at  some  specified  flow 
location.  This  condition  gives  rise  to  an  additional  equation  which  allows  the 
simultaneous  prediction  of  the  boundary  layer  development  and  the  streamwise 
static  pressure  distribution.  From  a straight-forward  applicaticn  of  the  con- 
tinuity equation,  it  can  be  shown  that  within  the  boundary  layer  .he  flow 
angle,  0,  at  a distance  T]  from  the  wall  is  given  by 


where  6*  the  displacement  thickness,  p the  edge  (wake)  density,  the  edge 
(wake)  velocity,  a the  wall  angle,  and  x is  the  principal  flow  direction.  Using 
Eq.  (l5)  the  inner  and  outer  flows  can  be  matched  anywhere  between  the  wall  and 
the  boundary  layer  edge;  however,  by  matching  at  the  displacement  surface 
(4  = 6*),  the  second  term  on  the  right-hand  side  of  Eq.  (15)  becomes  identically 
zero  and  the  resulting  expression  becomes 


tan  9 


r*  - 


As  can  be  seen  from  Eq.  (15),  the  choice  of  location  at  which  the  flows 
are  to  be  matched  is  somfewhat  arbitrary.  Two  locations  which  have  been  used 
commonly  are  the  boundary  layer  edge,  6,  ana  the  displacement  surface,  6*. 

For  relatively  thick  boundary  layers,  either  choice  represents  a compromise. 
If  6 matching  is  chosen,  the  entire  viscous  region  is  included  below  the 
matching  line;  however,  the  region  below  the  matching  line  also  includes  a 
region  of  transverse  pressure  gradients  which  are  ignored  in  the  boundary 
layer  calculation.  When  6*-matching  is  used,  the  region  below  the  matching 
line  is  expected  to  contain  negligible  transverse  pressure  gradients;  however, 
significant  viscous  and  rotational  effects  may  be  present  above  the  matching 
line.  In  the  present  effort  6*-matching  was  used. 
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The  boundary  layer  equations  are  solved  by  a forward  marching  procedure 
in  which  the  inner  and  outer  flows  are  coupled  at  each  step.  The  procedure 
is  an  iterative  one  which  first  approximates  the  displacement  thickness  (or 
edge  flow  angle  from  Eq.  (l6))  through  an  extrapolation  of  the  solution  a4 
previous  stations.  The  boundary  layer  momentum  equation  then  is  solved  using 
the  specified  displacement  thickness  to  obtain  both  the  boundary  layer  flow 
field  and  the  edge  flow  conditions.  Thus,  for  an  assumed  flow  angle,  + he 
viscous  calculation  predicts  a certain  edge  velocity.  Using  an  outer  inviscid 
flow  relation  and  the  boundary  layer  edge  velocity,  this  flow  angle  is  deter- 
mined again  and  the  two  values  of  the  flow  angle  are  compared.  If  the  two 
flow  angles  agree  within  a tolerance,  the  outer  and  inner  flows  are  considered 
to.be  matched.  If  not,  a new  displacement  thickness  is  assumed  and  the  proce- 
d’vire  repeated.  'The  modifications  in  the  equations  which  allow  specification 
of  6*  for  the  boundary  layer  solution  are  presented  in  Appendix  B. 

After  the  momentum  equation  is  matched  with  the  outer  flow  law,  the  energy 
equation  is  solved.  The  velocity,  density,  and  temperature  field  are  then 
evaluated.  The  matching  procedure  is  repeated  at  each  streamwise  step  and 
thus,  a solution  is  produced  by  marching  in  the  streamwise  direction. 

Outer  Flow  Analysis 

As  discussed  above,  the  strong  interaction  procedure  requires  an  inviscid 
flow  law  which  gives  a unique  relation  between  the  flow  angles  and  velocity. 

In  the  present  study,  the  well-known  Prandtl-Meyer  relationship  for  flows 
with  waves  of  one  family  was  used  to  describe  the  outer  flow.  Under  the 
Prandtl-Meyer  relationship,  the  flow  angle,  u),  and  the  Mach  number,  M,  are 
related  by 


uj  - y/-y  | (M2  - I)  - tan-' v^M2  - I (17) 

where  y is  the  ratio  of  specific  heats.  At  hypersonic  Mach  numbers  the  tangent- 
wedge  relation  can  readily  be  substituted  for  the  Prandtl-Meyer  relation  and, 
in  fact,  if  desired,  the  viscous  procedure  can  be  coupled  to  a supersonic 
inviscid  flow  code  to  predict  as  much  of  the  outer  flow  field  as  desired. 

Separated  Flow 


Since  the  boundary  layer  equations  are  solved  by  a forward  marching 
procedure,  numerical  instabilities  are  expected  to  be  generated  in  the  reverse 
flow  region  when  the  reverse  velocities  are  large  or  when  the  reverse  flow 
region ^xtends  over  a .large  portion  of  the  flow.  As  shown  by  several  authors 
(e.g.,  h*f.  'll),  the  so\irce  of  this  instability  is  in  the  streamwise  convec- 
tive term.-i.  In  the  present  effort,  the  instability  is  suppressed  by  the  usual 
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method  of  approximating  the  convective  terms  in  the  reversed  flow  region.  In 
so  far  as  the  energy  equation  is  concerned,  the  streamwise  convective  term 
was  simply  neglected  in  regions  of  reversed  flow.  Since  the  present  solution 
of  the  boundary  layer  equations  is  based  upon  a stream  function  formulation, 
the  suppression  of  the  convective  instabilities  in  the  momentum  equation  may 
not  be  obvious  and,  therefore,  the  suppression  of  this  term  in  the  momentum 
equation  is  now  discussed  in  some  detail. 

In  the  present  procedure  a stream  function  formulation  is  used  in  which 
the  stream  function,  F,  is  of  the  form: 

dF  , P u 

dy  ' ' Peue  (l8) 

and  the  transverse  velocity  v is  related  to  the  stream  function  through  applica- 
tion of  the  continuity  equation 


dp  u dp\i 

d x + dy  :° 

which  yields  a relation  of  the  form: 
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In  Eqs.  (l8)  through  (20)  x and  y are  the  streamwise  and  transverse  coordinates, 
u and  v are  the  velocity  components  in  the  x and  y directions,  o is  the  density, 
and  the  subscript  'e1  indicates  the  quantity  is  evaluated  at  the  edge  of  the 
boundary  layer.  After  a coordinate  transformation  of  the  form 


x'  = x t)  - y I h*  (21) 

is  made,  the  partial  differential  equations  are  reduced  to  ordinary  differ- 
ential equations  by  approximating  derivatives  in  the  streamwise  direction  by 
two  point  backward  differences;  i.e.,  for  the  variable  F, 

dF  ( X , -iy  ) F(X  ,7f)- F(x-  Ax,q) 

dx  = Ax  (22) 

As  can  be  seen  from  the  original  form  of  the  momentum  equation 


du  dp  + dr 
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and  an  examination  of  Eq.  (20)  derivatives  of  the  stream  function,  F,  with 
respect  to  the  streamwise  coordinate  x arise  both  through  the  terms  pudu/dx, 
which  controls  streamwise  convection  andpv6u/6x,  which  controls  transverse 
convection.  When  contributions  of  both  these  toms  to  SF/'&x  were  neglected 
in  regions  of  reversed  flow,  the  solution  became  unstable.  However,  when 
only  the  contribution  of  pudu/ dx  was  neglected  in  the  reversed  flow  region, 
a stable  solution  was  obtained.  In  this  latter  case  the  factor  &F/&x  which 
arises  from  the  term  pvdu/ify  is  treated  as  a nonlinear  coefficient  whose 
value  was  approximated  by  extrapolation  of  the  stream  function  from  previous 
stations.  This  method  of  approximating  the  convective  terms  in  the  momentum 
and  energy  equations  rendered  the  system  marginally  stable  in  the  presence 
of  reversed  flow.  However,  when  calculating  turbulent  boundary  layers  with 
large,  high  speed  reversed  flow  regions  instabilities  occurred.  Following 
Reyhner  and  Flugge-Lotz  (Ref.  4l),  rather  than  simply  neglecting  these  con- 
tributions to  the  convective  terms  in  the  momentum  equations,  the  absolute 
value  of  the  appropriate  coefficients  were  used  to  provide  added  stability. 

Turbulence  Model 

In  the  early  stages  of  the  present  work  the  turbulence  model  used  was 
identical  to  that  used  in  the  weak  interaction  solution  (See  Appendix  C). 

Since  this  turbulence  model  was  developed  for  attached  flows,  the  resulting 
poor  agreement  for  separating  flows  between  calculations  and  data  was  not 
unexpected.  The  weak  interaction  turbulence  model  was  based  upon  a mixing 
length  model  in  which  the  mixing  length  was  determined  either  by  an  equi- 
librium turbulence  assumption  or  by  the  solution  of  the  integral  turbulence 
kinetic  energy  equation.  This  turbulence  model  was  well-known  and  when  used 
in  the  turbulence  kinetic  energy  mode  has  been  proven  accurate  for  a broad 
spectrum  of  attached  boundary  layer  flows  (Refs.  25,  33.  34  and  46).  In 
strong  interaction  separating  flows,  however,  the  basic  model  developed  for 
attached  flow  yielded  poor  agreement  with  data  even  when  used  in  the 
turbulence  energy  mode  and,  in  fact,  calculations  made  with  the  original  model 
in  the  turbulence  energy  mode  showed  little  difference  from  calculations  made 
with  the  original  equilibrium  model.  Thus,  both  models  appeared  inadequate 
in  separated  flows. 


In  a recent  pap^r  (Ref.  4o)  Owen  showed  that  in  the  case  of  confined  coaxial 
,ie*  s with  recircula*  ion,  the  turbulence  fol  lowed  the  dividing  streamline  and  appeared 
to  H ffuse  from  this  str^aml  ine  as  the  main  flow  progressed  downstream.  Little  turbu- 
lence was  found  in  t he  "eversed  flow  region.  Rased  upon  some  computational  turbulence 
s turtles  and  the  laf  a of  Over.,  a nodi  float  ion  of  the  basic  turbulence  model  was  postu- 
lated for  boundary  layers  containing,  recirculating  flow.  The  details  of  the  modi- 
fied turbulence  model  are  presented  in  the  Results  section  of  this  report  along 
with  comparison:'  of  computations  using  this  modified  model  and  test  data. 
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Solution  Of  The  Strong  Interaction  Boundary  Layer  Equations 


The  strong  interaction  boundary  layer  procedure  was  formed  hv  modifying 
the  existing  weak  interaction  boundary  layer  code  in  the  manner  described 
previously.  However,  the  application  of  the  forward  marching  strong  inter- 
action procedure  to  compression  corner  or  incident  shock  problems  requires 
an  iterative  method  of  solution  based  upon  the  branching  solutions  obtained 
when  the  interaction  is  treated  as  a forward  marching  problem.  Therefore, 
the  present  section  discusses  characteristics  of  the  family  of  branching 
solutions  that  can  emerge  in  strong  interaction  computations.  As  will  be 
discussed  in  detail  subsequently,  the  forward  marching  procedure  can  be 
posed  so  that  the  problem  contains  a free  parameter.  Different  choices  of 
the  free  parameter  lead  to  different  branching  solutions  and  the  iterative 
solution  consists  of  determining  what  value  of  the  free  parameter  is 
required  to  satisfy  specified  downstream  boundary  conditions.  Since  proper 
choice  of  the  free  parameter  produces  the  desired  computational  solution, 
characteristics  of  the  selection  and  iterative  updating  of  the  free  parameter 
also  are  discussed. 

Characteristics  of  the  Branches 


The  nature  of  the  solution  to  the  strong  interaction  boundary  layer 
equationshas  been  the  subject  of  a series  of  investigations  (Refs.  5,  19,  22, 
24,  27  and  4l).  In  brief  these  investigations  have  shown  that  the  strong 
interaction  boundary  layer  equations  can  be  satisfied  by  a family  of  solutions 
(see  Fig.  l)  consisting  of  the  classical  weak  interaction  solution  (termed 
the  fundamental  solution)  and  in  addition,  a family  of  branching  or  departure 
solutions  (termed  free  interaction  solutions).  From  a numerical  point  of 
view  Tyson  (Ref.  53)  has  shown  that  these  branching  solutions  can  emerge  in 
solving  finite  difference  equations  when  the  streamwise  step  size  is  chosen 
to  be  less  than  a characteristic  departure  length;  for  large  streamwise  step 
sizes  only  the  fundamental  weak  interaction  solution  emerges.  The  validity 
of  these  branches  as  eigensolutions  of  the  laminar  boundary  layer  equations 
was  demonstrated  by  Hankey,  Dwoyer,  and  Werle  (Ref.  20)  and  the  role  of  the 
initial  conditions  in  determining  which  branch  of  the  family  of  solutions 
emerges  was  investigated  by  Werle,  Dwoyer,  and  Hankey  (Ref.  54).  Both  Refs. 

20  and  54  demonstrate  that  a single  free  parameter  exists  in  the  initial 
conditions  to  the  solution  of  the  laminar  strong  interaction  boundary  layer 
equations  that  controls  which  of  the  many  physically  possible  branching 
solutions  will  emerge.  By  controlling  the  free  parameter  one  can  select  the 
branch  to  emerge  from  the  solution.  More  recently  Bertke,  Werle,  and  Polak 
(Ref. 5 ) have  demonstrated  that  the  strong  interaction  turbulent  boundary 
.Layer  equations  contain  branching  characteristics  analogous  to  the  laminar 
branches  described  above. 
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Downstream  Boundary  Conditions 


An  obvious  problem  which  arises  when  the  interaction  flow  field  is 
solved  via  a forward  marching  procedure  is  the  specification  of  +he  down- 
stream boundary  condition.  At  the  upstream  boundary  of  the  calculation 
region,  the  boundary  layer  velocity  and  density  distributions  are  specified. 

At  the  downstream  station  the  correct  boundary  condition  must  be  based  upon 
the  pressure  (or  Mach  number)  being  consistent  with  the  incident  shock 
strength  or  compression  corner  angle.  In  the  present  effort  the  condition 
chosen  is  that  the  axial  pressure  gradient  be  zero  when  the  static  pressure 
is  equal  to  the  inviscid  downstream  value. 

Since  the  strong  interaction  boundary  layer  equations  are  solved  by  a 
forward  marching  technique,  an  iterative  solution  procedure  is  employed  to 
couple  the  downstream  boundary  conditions  into  the  solution.  First,  values 
of  the  boundary  layer  properties  as  well  as  the  free  parameter  controlling 
the  branching  are  set  at  the  upstream  boundary  of  the  flow  field.  The 
boundary  layer  equations  coupled  to  an  outer  flow  law  are  then  integrated 
downstream  by  the  implicit  finite  difference  procedure  outlined  previously 
and  the  pressure  distribution  generated  by  the  solution  is  monitored.  If 
the  axial  pressure  gradient  is  zero  at  the  station  where  the  static  pressure 
is  that  given  by  inviscid  considerations  downstream  of  the  shock,  then  the 
problem  is  solved;  if  not,  the  free  parameter  at  the  initial  station  is 
changed  and  the  boundary  layer  equations  integrated  again.  The  iteration  on 
the  free  parameter  at  the  initial  value  plane  is  continued  until  the  down- 
stream boundary  conditions  are  met  to  the  desired  tolerance.  This  procedure 
is  simply  a varient  of  the  well  known  shooting  technique  for  two-point 
boundary  value  problems. 

Branch  Control 


The  boundary  layer  calculations  were  started  using  the  weak  interaction 
boundary  layer  procedure  to  predict  the  development  of  a boundary  layer  under 
an  imposed  constant  static  pressure  distribution.  After  marching  several 
steps  to  settle  out  the  effect  of  intial  conditions,  the  weak  interaction 
calculation  was  replaced  by  a strong  interaction  calculation.  At  the  point 
of  switching  to  the  strong  interaction  calculation  the  solution  branch  was 
selected  by  imposing  a pressure  gradient  at  the  last  weak  interaction  station. 

In  general,  the  stronger  the  adverse  pressure  gradient  imposed,  the  sooner 
the  compressive  branch  emerged.  The  more  favorable  the  imposed  gradient,  the 
later  the  compressive  branch  emerged.  Kven  more  favorable  imposed  pressure 
gradients  produce  expansive  branches.  Since  the  branching  behavior  is 
determined  by  the  details  of  the  boundary  layer  profiles  at  the  start  of  the 
interaction  (Ref.  5M,  it  was  not  possible  to  determine  in  advance  the  imposed 
pressure  gradient  that  would  produce  a given  branch.  Tn  general,  a series  of 
runs  were  made  for  each  case  to  determine  empirically  the  details  of  the 
relation  between  the  imposed  pressure  gradient  and  the  emerging  branch  solution. 
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The  free  parameter  that  was  used  to  impose  the  pressure  gradient  at 
the  start  of  the  strong  interaction  was  a perturbation  in  the  edge  velocity 
termed  a "kick”.  Obviously,  the  perturbation  in  edge  velocity  imposed  a 
specified  pressure  gradient  upon  the  boundary  layer  at  the  last  weak  inter- 
action station.  Although  an  imposed  pressure  gradient  was  used  as  the  "free 
parameter"  in  the  present  study,  this  was  only  one  of  several  possible 
methods  which  could  be  used  to  generate  the  required  series  of  branching 
solutions.  In  the  process  of  iterating  on  a solution,  small  changes  in 
the  kick  at  the  start  of  the  interaction  are  required  to  produce  changes 
in  the  flow  field  at  the  downstream  boundary.  However,  in  the  process  of 
integrating  the  strong  interaction  boundary  layer  equations  downstream, 
small  changes  in  the  initial  flow  field  become  smothered  by  truncation 
error  from  the  linearizations  and  from  the  classical  problem  of  significant 
figures.  One  solution  to  this  problem  has  been  to  go  to  double  precision 
on  the  computer  and  thereby  gain  additional  significant  figures  in  the 
calculation . 

Another  method  for  controlling  the  branches  was  found  that  obviated 
the  need  for  double  precision.  This  method  consists  of  dividing  the 
interaction  procedure  into  two  parts,  each  using  a different  but  conceptually 
equivalent  free  parameter.  First  the  kick  is  used  as  the  free  parameter 
that  determines  the  location  of  the  start  of  the  free  interaction  relative 
to  the  location  of  the  incident  shock  or  compression  corner.  The  iteration 
on  the  kick,  described  above,  is  used  until  the  limit  in  significant  figures 
was  approached.  The  kick  is  then  frozen  at  the  latest  value.  The  second 
part  of  the  iteration  uses  the  distance  between  the  start  of  the  free 
interaction  and  the  location  of  the  incident  shock  or  compression  corner 
as  the  free  parameter.  Leaving  the  upstream  solution  unchanged,  the  location 
of  the  shock  of  corner  is  moved  up  or  downstream  until  the  downstream  flow 
field  is  compatible  with  the  downstream  boundary  conditions.  The  change  in 
free  parameter  used  to  control  the  iteration  has  removed  the  need  for  double 
precision  in  the  cases  that  were  investigated. 

Results 

The  strong  interaction  boundary  layer  was  computed  by  iteratively 
selecting  compressive  branches  of  the  boundary  layer  equations  until  compat- 
ibility with  the  downstream  boundary  conditions  was  obtained.  This  section 
presents  the  results  of  applying  this  strong  interaction  boundary  layer 
calculation  procedure  to  a series  of  laminar  and  turbulent  compression  corner 
flows.  The  laminar  runs  are  presented  to  verify  the  procedure  used  since 
the  laminar  calculations  are  not  complicated  by  turbulence  modeling.  Both 
the  validity  of  the  numerics  and  the  control  of  the  branching  solutions  are 
presented.  When  turbulent  strong  interaction  calculations  were  performed, 
the  qualitative  behavior  of  the  solutions  was  found  to  be  the  same  as  in  the 
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laminar  case,  however  originally  the  quantitative  agreement  with  experimental 
data  was  found  to  be  poor.  The  source  of  the  problem  in  the  case  of  turbulent 
flows  was  found  to  be  the  turbulence  model.  Therefore,  based  upon  some  recent 
experimental  data  of  Owen  (Ref.  40  ) a revised  turbulence  model  was  postulated. 
Turbulent  solutions  with  both  the  original  and  the  revised  model  are  presented 
subsequently. 

Laminar  Flow 


Other  investigators  have  achieved  apparently  accurate  predictions  of 
laminar  interacting  boundary  layers  using  physical  assumptions  similar  to 
those  used  in  the  present  effort.  These  available  solutions  can  be  used  to 
confirm  the  correct  operation  of  the  UTRC  Code.  Furthermore,  the  laminar 
case  provides  a direct  confirmation  of  the  numerics  since  it  is  not  complicated 
by  the  turbulence  modeling  required  for  turbulent  boundary  layer  calculations. 

In  separating  boundary  layers,  experimental  data  (Refs.  11,  13,  l6  and  4-9) 
indicates  that  when  separation  in  supersonic  flow  results  from  an  incident 
shock  or  a compression  corner,  the  shape  of  the  initial  pressure  rise  and 
the  associated  changes  in  the  boundary  layer  profile  are  independent  of  the 
specific  cause  of  the  compression,  a phenomenon  termed  free  interaction. 
Experimental  data  also  shows  that  the  distance  between  the  start  of  the 
pressure  rise  and  the  incident  shock  or  compression  corner  increases  as  the 
strength  of  the  imposed  compression  increases.  Thus  for  given  conditions, 
free  interaction  pressure  distributions  form  a family  of  identical  initial 
pressure  rise  contours  (Fig.  2)  whose  location  relative  to  the  location  of 
the  incident  shock  or  compression  corner  is  determined  bv  the  strength  of 
the  overall  compression.  Two  free  interaction  pressure  distributions  are 
parallel  to  each  other  only  being  separated  by  a constant  displacement  in 
the  streamwise  direction.  It  should  be  noted  that  both  laminar  and  turbulent 
boundary  layers  have  been  found  experimentally  to  exhibit  this  free  inter- 
action property.  In  the  case  of  laminar  flow  the  various  pressure  rise 
curves  reach  a common  plateau  at  a pressure  level  termed  the  "plateau  pressure" 
(see  Fig.  2)  which  occurs  at  a streamwise  location  upstream  of  the  incident 
shock  or  compression  corner. 

In  the  present  effort  in  which  strong  interacting  boundary  layers  ere 
computed  by  a forward  marching  procedure,  departure  solutions  emerge  and,  as 
described  in  the  previous  section,  these  computational  branching  solutions  are 
bona  fide  solutions  of  the  boundary  layer  equations.  The  behavior  of  these 
computational  branching  solutions  was  found  to  correspond  in  a qualitative 
manner  to  the  experimentally  observed  free  interactions  for  both  laminar  and 
turbulent  separating  boundary  layers.  For  given  flow  properties  different 
departure  solutions  give  pressure  rises  which  initially  are  parallel  and  only 
differ  in  the  location  of  the  pressure  rise.  The  computed  laminar  branching 
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Figure  2.  Experimental  Free  Interaction  Pressure  Rise 


o*.  mu  i 


solutions  obtained  in  the  present  study  were  compared  to  the  free 
interacting  plateau  pressure  correlation  of  Curie  (Ref.  14)  and  the  pressure 
rise  to  separation  correlation  of  Chapman,  et  al  (Ref.  11) in  Fig.  3.  The 
comparisons  of  the  calculations  with  the  data  correlations  show  the 
predictions  of  the  present  procedure  to  be  in  good  agreement  with  both  the 
separation  pressure  and  the  pleateau  pressure  correlations. 

Two  calculations  corresponding  to  two  different  upstream  conditions 
are  shown  in  Fig.  3*  In  both  calculations  the  interacting  solutions  are 
started  by  first  marching  the  code  for  several  streamwise  stations  under 
a specified  free-stream  velocity  distribution,  i.e.,  an  imposed  pressure 
distribution.  A velocity  perturbation,  the  "kick", is  then  added  to  the 
imposed  edge  velocity  to  trigger  the  branching  solutions.  At  stations  down- 
stream of  the  kick,  the  edge  velocity  is  determined  from  the  strong- 
interaction  calculation . Different  kicks  lead  to  different  profiles  at 
the  last  station  prior  to  the  interaction  and  thus  lead  to  different- 
solution  branches.  The  desired  branch  can  be  selected  by  iteratively  changing 
the  "kick".  This  control  of  the  branches,  along  with  the  ability  to  calculate 
reattachment  was  confirmed  by  running  a series  of  sample  calculations  at 
M = 4.  The  branch  control  obtained  via  an  upstream  kick  is  presented  in 
Fig.  4. 

When  the  conditions  at  the  start  of  the  interaction  are  varied  in 
an  attempt  to  match  downstream  boundary  conditions,  small  changes  in  the 
kick  are  required.  Because  of  the  tolerances  used  in  the  iteration  and  the 
truncation  inherent  in  the  solution,  the  effect  of  small  changes  in  the  "kick", 
smaller  than  the  changes  shown  in  Fig.  4,  can  become  numerically  insignificant 
by  the  time  the  calculation  marches  past  the  corner.  Tighter  tolerances  and 
double  precision  could  be  used  to  alleviate  this  problem,  however  either  of 
these  would  have  an  associated  increase  in  computer  run  times.  Therefore,  a 
second  method  of  branch  control  was  investigated  in  order  to  control  the 
reattachment  and  downstream  branching.  The  second  method  of  control  was 
obtained  by  displacing  the  compression  corner  small  distances  up  or  down- 
stream relative  to  the  start  of  the  interaction.  Control  of  the  reattachment 
and  downstream  branching  solutions  was  successfully  obtained  by  this  means  as 
shown  in  Fig.  5.  In  practice  a combination  of  both  methods  of  branch  control 
is  employed.  The  "kick"  in  edge  velocity  at  the  start  of  the  intex-action  is 
used  to  control  the  branching  until  the  downstream  solution  became  insensi- 
tive to  small  changes  in  the  "kick".  Final  resolution  of  the  downstream 
boundary  condition  is  accomplished  by  small  displacements  of  the  location  of 
the  compression  corner.  Applying  this  two  stage  method  of  branch  control, 
a laminar  strong  interaction  calculation  was  made  to  compare  to  the  data  of 
Lewis,  Kubota,  and  Lees  (Fig.  6).  After  a suitable  kick  was  imposed,  various 
branching  solutions  were  obtained  by  corner  displacement.  As  expected,  two 


Figure  3.  Laminar  Free  Interaction  Calculations 


families  of  solutions  emerged:  one  giving  continually  increasing  pressure 

and  one  giving  a pressure  maximum  followed  by  a pressure  decrease.  Since 
this  study  is  primarily  concerned  with  strong  interacting  turbulent  boundary 
layers,  the  pressure  overshoot  that  occurs  in  laminar  boundary  layers  after 
reattachment  (Ref.  28)  is  not  addressed  in  this  report.  The  corner  location 
was  iterated  upon  to  move  the  branching  location  between  the  two  families 
further  and  further  downstream.  The  results  of  the  calculation  are  presented 
in  Fig.  6 which  shows  good  agreement  with  data  over  most  of  the  interesting 
range. 

Turbulent  Flow 


After  the  procedure  was  confirmed  through  the  laminar  calculations 
presented  in  Figs.  3-6  the  procedure  was  applied  to  strong  interacting 
turbulent  boundary  layers.  The  results  appeared  to  be  qualitatively  reasonable, 
however,  the  free  interaction  pressure  rise  was  more  abrupt  and  led  to 
higher  plateau  pressure  than  was  indicated  by  data  as  shown  in  Fig.  7. 

Similar  results  were  obtained  by  Bertke,  Werle,  and  Polak  (Ref . 5 ) using 
a different  turbulence  model.  In  the  present  code  the  computed  turbulence 
structure  is  based  upon  a mixing  length  model  in  which  the  mixing  length 
can  be  determined  by  either  an  equilibrium  turbulence  assumption  (equilibrium 
model)  of  by  the  solution  of  the  integral  turbulence  kinetic  energy  equation 
(TKE  model).  Under  either  the  equilibrium  or  TKE  option,  the  turbulence 
structure  is  based  upon  an  assumed  one-parameter  mixing  length  profile  which 
varies  with  distance  from  the  wall  as  a hyperbolic  tangent  function.  In 

the  immediate  vicinity  of  the  wall  the  mixing  length  is  damped  by  a prob-  , ; 

ability  damping  function.  The  free  parameter  in  the  profile  is  the  "wake" 

value  of  the  mixing  length;  i.e.,  the  value  far  removed  from  the  wall.  When 

the  equilibrium  turbulence  option  is  used,  this  wake  value  is  a function  of 

the  boundary  layer  thickness  and  the  local  momentum  thickness  Reynolds  number. 

When  the  turbulence  kinetic  energy  option  is  used,  the  "wake"  value  of  the 
mixing  length  emerges  from  the  solution  of  the  integral  turbulence  kinetic 
energy  equation. 

A turbulent  free  interaction  was  computed  first  with  the  equilibrium 
model  and  then  with  the  TKE  model.  As  shown  in  Fig.  8,  the  calculated  wall 
pressure  and  skin  friction  distributions  were  found  to  be  virtually  identical 
even  though  the  wake  mixing  length  distributions  of  the  two  methods  differed 
considerably  (see  Fig.  9).  The  low  Reynolds  number  correction  factor 
(Appendix  C)  is  used  in  both  cases.  In  an  attempt  to  gain  some  insight  into 
the  sensitivity  of  the  wall  pressure  distribution  to  the  turbulence  modeling 
assumptions  a series  of  frozen  turbulence  models  were  evaluated.  First,  the 
eddy  velocity  was  frozen  at  the  value  at  the  start  of  the  interaction  and 
held  constant  along  lines  of  constant  distance  from  the  wall  (y  = const). 

Second,  the  turbulence  intensity  (Reynolds  shear  stress  -u'v')  was  frozen 
at  the  start  of  the  interaction  and  again  held  constant  along  y = const. 
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Third,  the  eddy  viscosity  was  frozen  and  held  constant  along  streamlines, 
and  fourth,  the  turbulence  intensity  was  frozen  and  held  constant  along 
streamlines . 

The  results  presented  in  Fig.  10  demonstrate  that  the  wall  pressure 
distributions  could  be  drastically  altered  by  the  turbulence  modeling 
assumptions.*  However,  as  shown  in  Figs.  9 and  10  the  wall  pressure  distri- 
butions produced  by  the  original  TKE  and  equilibrium  models  were  very  similar 
yet  the  wake  mixing  length  distributions  produced  by  these  assumptions  were 
considerably  different.  To  resolve  this  apparent  contradiction  an  inves- 
tigation was  made  into  the  details  of  the  turbulence  structure  produced  by 
the  TKE  and  equilibrium  models.  The  mixing  length  distributions  within  the 
boundary  layer  produced  by  the  TKE  and  equilibrium  models  were  compared  as 
shown  in  Fig.  11  at  three  streamwise  locations.  Although  the  TKE  and 
equilibrium  models  produce  significantly  different  values  of  the  wake  mixing 
length,  Fig.  11  shows  that  the  turbulence  modeling  near  the  wall  remains 
unaffected  by  the  changes  in  the  wake  turbulence.  Thus  the  problem  with  the 
turbulence  modeling  may  not  be  one  of  predicting  the  wake  properties  but  may 
be  one  of  predicting  the  distribution  of  turbulence  near  the  wall,  especially 
in  the  recirculating  region. 

In  a recent  paper  (Ref.  kO  ) Owen  showed  that  in  the  case  of  confined 
coaxial  jets  with  recirculation,  the  turbulence  followed  the  dividing  stream- 
line separating  the  recirculating  flow  from  the  outer  flow  (in  a time  averaged 
sense)  and  appeared  to  diffuse  from  this  streamline  toward  the  wall  as  the 
main  flow  progressed  downstream.  Little  turbulence  was  found  in  the  reversed 
flow  region.  Based  upon  the  results  of  the  turbulence  studies  presented  in 
Figs.  9“H  and  upon  the  data  of  Owen,  a modification  of  the  basic  turbulence 
model  in  the  near  wall  region  has  been  postulated  for  boundary  layers  containing 
recirculating  flow.  This  modification  employs  the  observation  that  in  attached 
boundary  layers  the  wall  is  a continuous  streamline  of  the  flow.  In  separated 
boundary  layers,  the  dividing  streamline  is  the  limiting  streamline  of  the 
main  flow  region,  and  the  turbulence  approaching  the  recirculation  region  is 
convected  around  the  bubble  by  the  mean  flow.  Therefore,  the  damping  of  the 
turbulence  from  the  wall  in  the  original  turbulence  model  is  replaced  by 
damping  of  the  turbulence  from  the  dividing  streamline  when  a recirculation  is 
present.  Since  Owen's  data  implies  a slight  spreading  of  the  turbulence  into 
the  recirculating  region,  a diffusion  of  the  turbulence  from  the  dividing 
streamline  into  the  recirculating  flow  region  is  included  in  the  modified 
model.  This  diffusion  is  postulated  to  spread  the  turbulence  at  a U°  half  angle 
into  the  forward  flow  portion  of  the  recirculating  region.  Damping  by  the  wall 
is  expected  to  have  a strong  effect  in  the  relatively  low  speed  reversed  flow, 
and  since  Owen  found  little  turbulence  in  this  region  the  turbulence  is  assumed 
to  be  completely  damped  in  the  reversed  flow.  In  attached  boundary  layers,  the 
modified  turbulence  model  collapses  into  the  original  model.  Mikulla  and 
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Figure  9 . Wake  Mixing  Length  Distribution  With  Attached  Boundary  Layer 
Turbulence  Model  Including  Low  Reynolds  Number  Correction 
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THEORY 


1 EQUILIBRIUM  TURBULENCE 

2 TURBULENCE  KINETIC  ENERGY  EQUATION 

3 FROZEN  VISCOSITY  ALONG  Y = CONST. 

4 FROZEN  TURBULENCE  ALONG  Y CONST. 

5 FROZEN  VISCOSITY  ALONG  S.L. 

6 FROZEN  TURBULENCE  ALONG  S.L. 

DATA  OF  SPAID  AND  FRISHETT  (REF  51) 


Horstman  (Ref.  36)  have  recently  published  data  that  includes  turbulence 
measurements  for  a shock  wave-boundary  layer  interaction.  Although  the 
measurements  taken  only  roughly  define  the  turbulence  field,  the  data  is 
in  qualitative  agreement  with  the  conclusion  drawn  from  the  data  of  Ref.  40. 

The  turbulent  free  interaction  calculation  of  Fig.  7 was  recalculated 
using  the  modified  turbulence  model.  A comparison  of  the  wall  pressure 
distributions  calculated  using  the  original  and  the  modified  turbulence 
model  is  presented  in  Fig.  12,  where  the  wake  mixing  length  has  been  set 
in  both  cases  by  the  Prandtl  mixing  length  equilibrium  turbulence  model. 

As  shown  in  Fig.  12  the  modified  turbulence  model  corrects  the  overly  abrupt 
initial  pressure  rise  and  reduces  the  plateau  pressure,  and  therefore,  the 
modified  model  leads  to  better  agreement  with  experimental  data.  A comparison 
of  the  turbulence  distributions  for  the  original  and  modified  turbulence 
models  is  presented  in  Figs.  13-16.  It  can  be  seen  from  Figs.  13-16  that 
the  modifications  in  the  turbulence  models  affect  the  wall  region  and  in 
the  case  of  the  TKE  model  can  also  affect  the  wake  region. 

The  calculations  of  the  strong  interaction  boundary  layer  procedure 
with  the  modified  turbulence  model  are  compared  to  the  compression  corner 
data  of  fipaid  and  Frichett  (Ref.  51)  in  Figs.  17-19*  Three  cases  at  an 
upstream  Mach  number  of  2.93  were  investigated  for  ramp  angles  of  9*8l°» 

16.06°  and  I9.670.  The  separation  bubble  length  varies  in  these  cases  from 
very  small,  less  than  a third  of  the  upstream  boundary  layer  thickness,  to 
very  large,  over  2 l/2  times  the  upstream  boundary  layer  thickness.  Figures 
17-19  show  in  all  three  cases  that  the  calculated  overall  pressure  distribution 
is  in  good  agreement  with  the  data. 
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Figure  12.  Free  Interaction  Wall  Pressure  Distributions 
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Figure  13.  "Equilibrium"  Wake  Mixing  Length  Distributions  With  Low 
Reynolds  Number  Correction 
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Figure  18.  Turbulent  Strong  Interaction  Solution  For  16.06°  Ramp  Using  Modified  Turbulence  Model 
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TIE  NAVIER-STOKES  APPROACH 

The  second  phase  of  the  present  effort  is  a feasibility  study  in  which 
the  Multi-dimensional  Implicit  Nonlinear  Time -dependent  (MINT)  procedure  of 
Qriley  and  McDonald  (Ref.  9)  is  applied  to  the  interaction  flow  field.  Until 
fairly  recently  it  would  have  been  impractical  to  consider  a Navier-Stokes 
solution  to  this  problem,  however,  rapid  advances  both  in  numerical  analysis 
and  computer  technology  now  make  the  Navier-Stokes  procedure  a possible  alter- 
native to  the  extended  boundary  layer  procedures  for  shock  wave-boundary  layer 
interaction  predictions.  Solutions  based  upon  a Navier-Stokes  procedure  have 
several  advantages  over  those  based  upon  boundary  layer  analyses.  First  of 
all,  solutions  based  upon  the  Navier-Stokes  equations  need  not  make  an  arbi- 
trary division  between  viscous  and  inviscid  portions  of  the  interaction  flow 
field.  In  addition,  Navier-Stokes  procedures  solve  a full  transverse  momentum 
equation  and  need  not  make  any  approximation  to  the  convection  terms  in  sepa- 
rated flow  regions.  These  considerations  indicate  that  even  if  a Navier-Stokes 
interaction  procedure  would  require  more  computer  resources  (storage  and  run 
time)  than  a boundary  layer  interaction  procedure,  its  potentially  increased 
accuracy  may  make  it  an  attractive  alternative. 

In  the  past  few  years  several  investigators  have  applied  Navier-Stokes 
procedures  to  the  interaction  problem  (e.g. , Refs.  2,  29,  30,  47  and  48).  With 
the  exception  of  MacCormack  (Ref.  30),  all  those  solutions  have  been  based  upon 
explicit  finite  difference  procedures.  In  Ref.  30,  MacCormack  used  a hybrid 
procedure  in  which  convective  terms  were  treated  explicitly  and  viscous  terms 
treated  implicitly.  One  major  difficulty  which  seems  to  emerge  in  varying 
degrees  of  severity  in  Navier-Stokes  interaction  solutions  is  the  accurate 
treatment  of  the  shock  wave.  Over  the  past  years  a considerable  effort  has 
been  expended  in  developing  explicit  procedures  which  can  solve  the  flow 
equations  in  the  presence  of  flow  discontinuities  (shocks).  The  explicit 
finite  difference  interaction  procedures  mentioned  above  are  based  upon  such 
a procedure  originally  suggested  by  the  weak  solution  concepts  of  Lax  and 
developed  to  a fine  point  by  MacCormack;  this  procedure  is  described  in  detail 
in  Ref.  29.  In  a separate  development  Moretti  (Ref.  37)  treated  shock  waves  as 
spatial  discontinuities  with  explicit  recognition  of  their  presence.  However, 
at  the  present  time  this  technique  cannot  be  used  in  viscous  regions  where  the 
shock  can  be  physically  diffracted  by  the  shear  layer.  To  date,  no  similar 
large  scale  shock  treatment  effort  has  been  undertaken  using  implicit  methods. 

Tn  a study  of  inviscid  transonic  flow  Murman  and  Cole  (Ref.  38)  solved  the 
transonic  perturbation  equations  using  an  implicit  line  relaxation  procedure. 

The  Murman-Colo  technique  succeeded  in  capturing  shocks  in  supercritical  flow 
as  part  of  a continuous  solution  owing  largely  to  the  numerical  viscosity 
introduced  by  the  difference  scheme.  In  a later  work  Murman  (Ref.  39)  used  a 
special  difference  operator  at  the  sonic  line  (the  shock  point  operator')  and 
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again  captured  the  shock  as  part  of  a continuous  solution.  The  shock 
representation  in  this  latter  work  was  considerably  sharper  than  in  the  Murman 
and  Cole  work.  Both  procedures  are  discussed  in  some  detail  by  Hafez  and 
Cheng  (Ref.  19).  In  addition,  Beam  and  Warming  (Ref.  4)  have  used  the  Briley- 
McDonald  numerical  procedure  to  investigate  the  flow  field  about  a shock  in 
which  the  flow  passes  from  a supersonic  to  a subsonic  state.  Therefore, 
although  some  effort  has  been  expended  upon  shock  representation  within 
implicit  schemes,  the  effort  has  concentrated  primarily  upon  (i)  flows  which 
are  inviscid  and  (ii)  shocks  through  which  the  flow  goes  from  a supersonic  to 
a subsonic  state.  Neither  of  the  above  requirements  is  satisfied  in  the  shock 
wave-boundary  layer  interaction  problem  where  the  flow  is  viscous  and  the  shock 
is  usually  not  strong  enough  to  cause  subsonic  flow  on  its  downstream  side. 

The  Navier-Stokes  portion  of  the  present  effort  represents  a simple  feasibility 
study  in  which  a fully  implicit  procedure  is  applied  to  the  shock  wave-boundary 
layer  interaction  problem.  The  study  assesses  an  application  of  the  MINT  code, 
in  its  present  form,  to  the  interaction  problem,  gives  an  estimate  of  computer 
run  times  and  indicates  areas  in  which  further  development  work  would  be 
required  to  obtain  a reliable  and  accurate  prediction  procedure. 

The  Basic  Analysis 

The  Governing  Equations 


'The  equations  solved  in  the  procedure  represent  conservation  equations  of 
mass,  momentum  and  energy  (Ref.  6).  For  simplicity  the  equations  are  expressed 
in  vector  notation  below  and  all  quantities  are  nondimensional.  Velocities  are 
normalized  by  Up,  density  by  Pp,  enthalpy  by  hp,  temperature  by  Tp,  pressure  by 
Pp  = ppR^Tp  where  R is  the  gas  constant,  dynamic  viscosity  by  pp,  and  time  by 
(L/Up)  where  L is  tne  reference  length,  body  forces  and  bulk  viscosity  are  all 
assumed  to  be  negligible.  The  resulting  time-averaged  equations  are  given  by 


Continuity 


If  =-V(pu) 


(A) 


Conservation  of  Momentum 

d(PQ)  _ _ y (puiJ)  - ^ 

'd"d 


at 


= -?  ^ + ^el.5)  - f k*  h.f,V  0)]  < ’5' 


Conservation  of  Energy 

djp_ 
at 


_ V (puH)  + A-  + J_  v (f  VH) 


Poh0  at  Re 


(26) 


V(SJI)] 


42 


I 


fa. 


The  mean  flow  rate  of  strain  tensor  in  Eq.  (25)  is  given  by 

4 = -i-  [(Vo)+(Vu)T]  (27) 

The  necessary  thermodynamic  relationships  are 


P = PT 


H=  h + _L  Hei  (u  u) 
2 hp 

and  for  constant  specific  heat  the  enthalpy  is 


h = CpT 


(28) 

(29) 


(30) 


In  order  to  solve  the  above  system  of  equations  it  is  necessary  to  specify 
the  turbulent  exchange  coefficients  and  r^.  In  the  present  analysis  since 

the  effective  Prandtl  number  is  defined  from  knowledge  of  turbulent  flows  of 
gases,  only  the  turbulent  momentum  exchange  coefficient,  ucl'p.  must  be  speci- 
fied, The  energy  exchange  coefficient  is  obtained  from  the  relation 
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and  the  effective  viscosity  is  obtained  from  a turbulence  model. 

The  'Turbulence  Model 

In  the  case  of  laminar  flow  the  governing  equations,  Eqs.  (24)  - (26), 
along  with  the  relations  expressed  by  Eqs,  (27)  - (31)  are  sufficient  to  deter- 
mine a solution  when  proper  boundary  conditions  are  applied.  Ifowever,  in  tur- 
bulent flow  it  is  necessary  to  hypothesize  a turbulence  model  relating  the 
turbulent  viscosity  to  the  other  problem  variables.  Over  the  past  years  many 
such  models  have  been  hypothesized.  These  include  equilibrium  models  relating 
a mixing  length  or  eddy  viscosity  to  local  variables,  and  historical  models  in 
which  the  local  turbulent  stress  is  determined  through  an  additional  equation 
(or  equations)  relating  the  stress  to  the  upstream  history  of  the  flow.  A 
model  of  the  latter  type  in  which  a mixing  length  is  hypothesized  but  in  which 
the  magnitude  of  the  mixing  length  is  determined  by  the  turbulence  kinetic 
energy  equation  is  discussed  in  the  boundary  layer  section  of  this  report. 
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Since  the  present  Navier-Stokes  effort  represents  an  initial  feasibility- 
study,  a relatively  simple  turbulence  model  based  upon  a mixing  length  was 
used.  This  model  is  similar  to  the  two-region  eddy-viscosity  model  utilized 
by  Shang  and  Hankey  (Ref.  1*8),  in  which  the  mixing  length  in  the  inner  region 
is  given  by  the  equilibrium  model  discussed  in  Section  II,  i.e.,  the  turbulent 
viscosity  is 


My.  = pJL  2 V2£  e 


(3?) 


with  the  mixing  length,  l , given  by 
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where  6^  is  the  local  boundary  layer  thickness,  K is  the  von  Karman  constant, 
y is  the  distance  from  the  wall,  and  is  a sublayer  damping  factor  defined 
by 


1/2 


(33) 


where  P is  the  normal  probability  function,  and 

»**  y (f),/2  ir  m 

Here  t is  the  local  shear  stress,  y+  = 23  and  ct+  = 8. 

For  two-dimensional  equilibrium  turbulent  boundary  layers  has  a 

value  of  approximately  0.1,  which  was  employed  in  the  calculations  reported 
herein.  In  addition,  because  of  the  difficulty  of  defining  the  boundary  layer 
edge  in  an  interacting  flow  field,  6b  was  taken  as  the  boundary  layer  thickness 
upstream  of  the  interaction  region. 

The  outer  region  eddy  viscosity  is  given  by  the  Clauser  defect  law,  i.e., 

MTq  = 0.0168  />umax  8*nc  yk  (35) 
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where  6*pnc  is  the  kinematic  displacement  thickness 


■ 


(36) 


and  v is  the  Klebanoff  intermittency  which  may  be  approximated  by 
k 

H1  t55(fb)T'  <3' 

In  Eqs.  (35)  - (36),  umax  has  been  employed  instead  of  the  edge  velocity  to 
avoid  possible  anomalous  behavior  of  the  transient  solution  (Ref.  48). 

Artificial  Damping 


The  treatment  of  shock  waves  in  a Navier-Stokes  calculation  procedure 
with  finite  differences  requires  the  use  of  artificial  diffusion  in  order  to 
prevent  numerical  oscillation  which  may  cause  solution  divergence.  In  the 
present  implicit  procedure  a modified  form  of  the  fourth-order  pressure  damp- 
ing term  suggested  by  MacCormack  and  Baldwin  (Ref.  29)  has  been  employed, 
viz.,  a term  of  the  form. 
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(38) 


has  been  added  to  each  of  the  governing  equations  being  solved  for  each 
'•oordinate  direction  xk>  The  diffusion  coefficient  T)k  in  Eq.  (38)  was  taken 
as 
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where  0 < 8 < \ , c is  the  local  sound  speed,  and  the  average  pressure  p 
given  by 


P = 


+ 


(40) 


where  i is  the  grid  point  index  in  the  x^-direction.  In  the  damping  term  the 
factor  f is  set  to  1.6  for  the  continuity  equation  and  to  p for  all  other 
conservation  equations. 


Boundary  Conditions 

In  the  present  Navier-Stokes  calculations  boundary  conditions  were 
specified  along  all  boundaries  of  the  computational  regime.  At  the  wall  the 
normal  and  tangential  velocities  were  set  to  zero  and  a three  point  one-sided 
difference  form  of  the  continuity  equation  was  applied.  At  the  upstream  edge 
of  the  calculation  region  both  velocity  components  as  well  as  the  density  were 
specified.  These  values  were  determined  by  performing  a calculation  over  a 
region  extending  from  a station  slightly  downstream  of  the  plate  leading  edge 
to  a station  slightly  upstream  of  the  location  at  which  the  interaction  calcula- 
tion was  to  be  initiated.  This  calculation  yielded  a flat  plate  boundary  layer 
solution  that  was  consistent  with  the  finite  difference  representation,  turbu- 
lence model  and  grid  spacing  used  for  the  interaction  problem.  The  profiles 
calculated  at  the  station  having  the  desired  boundary  layer  thickness  then  were 
used  as  upstream  boundary  conditions  for  the  interaction  calculation.  At  the 
downstream  boundary  first  derivatives  of  all  dependent  variables  were  set  to 
zero.  Finally,  both  velocity  componentsand  the  density  were  set  at  the  outer 
edge  boundary.  This  boundary  was  taken  at  a distance  of  approximately  five 
boundary  layer  thicknesses  from  the  wall.  An  oblique  shock  was  assumed  to 
penetrate  the  boundary  between  the  third  and  fourth  grid  points  from  the 
upstream  station,  idle  values  downstream  of  the  shock  were  obtained  from  the 
oblique  shock  relations. 

The  Nonuniform  Grid 

The  accuracy  of  solutions  computed  with  a given  number  of  grid  points 
often  can  be  improved  by  using  a nonuniform  grid  spacing  to  ensure  that  grid 
points  are  closely  spaced  in  regions  where  the  solution  varies  rapidly.  In  the 
interaction  flow  field  large  gradients  are  present  near  the  wall  and  conse- 
quently fine  grid  resolution  is  desired  in  this  region.  This  grid  resolution 
was  obtained  using  an  analytic  coordinate  transformation  devised  by  Roberts 
(Ref.  4f)  which  is  a very  effective  means  of  introducing  a nonuniform  grid  when 
the  steep  gradients  occur  near  the  computational  boundaries.  Suppose  that  N 
gri  i points  are  to  be  used  in  the  range  X-^  £ X £ X-*.  and  that  steep  graidents 
are  anticipated  in  a region  of  thickness,  0 (X-,-Xj_)  near  X. . Then  Roberts' 
transformation  X..,(X)  is  given  by 

XT(X).N*(N-I)  (1*1) 
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where  a = X2  - X^,  b"  = a"/(l-8),  and  c » X,v  The  use  of  equally-spaced  points 
in  the  transformed  coordinate,  Xq>,  ensures  an  adequate  resolution  of  both  the 
overall  region  X-j_  £ X s X?  and  the  subregion  X1  s X s ^(Xj  -X1).  Derivatives 
with  respect  to  the  physical  coordinate,  X,  are  obtained  from  the  following 
formulas : 
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The  use  of  tliree-point  difference  operators  for  X.j,  derivatives  in  Eqs.  (33)  and 
(34)  produces  similar  operators  for  X derivatives.  These  X-derivative  operators 
can  be  computed  at  the  start  of  a calculation  and  stored,  along  with  the  X 
locations  of  grid  points. 

In  the  present  offer"  rather  than  use  the  Roberts  transformation  over  the 
entire  flow  regime,  the  transformation  was  used  only  in  the  vicinity  of  the 
wall.  In  regions  away  from  the  wall  an  equally  spaced  grid  was  used.  The  two 
regions  were  matched  by  requiring  that  at  the  join  point  the  physical  step  rise 
in  the  transformed  regime  be  the  same  as  the  step  size  in  the  equally  paced 
regime.  Thus  continuity  of  the  spacing  was  preserved.  It  should  be  noted  that 
although  the  mesh  was  continuous,  its  derivative  was  not  and  an  improved  mesh 
would  join  the  Roberts  region  to  the  equally  spaced  region  via  a transition 
region.  However,  the  mesh  did  not  seem  to  be  the  cause  of  any  numerical  prob- 
lems and,  therefore,  the  two  region  mesh  was  deemed  adequate  for  the  present 
investigation. 

Method  of  Solution 

Exact  analytical  solutions  of  the  Navier-Stokes  equations  are  rare  due  to 
their  high  order  and  coupled  nonlinearity.  As  a result  numerical  methods  must 
generally  be  employed  for  the  solution  of  these  equations.  In  addition,  one  of 
the  major  obstacles  to  the  solution  of  the  three-  imensional  compressible 
Navier-Stokes  equations  is  the  large  amount  of  computer  t ime  required,  and  con- 
sequently efficient  computational  methods  are  hiohl;.  miracle.  Most  numerical 
procedures  \ised  to  solve  the  Navier-Stokes  eqtiations  :;av  mvn  base  i on  explicit 
difference  schemes  for  the  unsteady  form  of  the  govrnii,.-  qua- ions,  and  are 
subject  to  one  or  more  stability  restrictions  on  the  size  of  the  time  step 
relative  to  the  spatial  mesh  size  (e.g..  Refs.  1,  '9  and  43).  These  stability 
limits  usually  correspond  to  the  well-known  Couran? -Trio  irichs-Lewy  (CFL)  con- 
dition and  in  some  methods  to  an  additional  viscous  stability  condition.  A 
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key  disadvantage  of  such  conditionally  stable  methods  is  that  the  maximum  time 
step  is  fixed  by  the  spatial  mesh  size  rather  than  the  physical  time  dependence. 
If  a steady  solution  is  being  computed  as  the  asymptotic  limit  of  the  unsteady 
solution, then  using  a small  time  step  requires  a large  number  of  steps  to  reach 
the  steady  solution. 

In  contrast  to  most  explicit  methods,  many  implicit  methods  tend  to  be 
stable  for  large  time  steps,  and  hence,  offer  the  prospect  of  substantial 
increases  in  computational  efficiency,  provided  of  course  that  the  computa- 
tional effort  per  time  step  is  competitive  with  that  of  the  conditionally 
stable  methods.  An  accurate  and  efficient  implicit  method  termed  the  MINT  pro- 
cedure has  been  developed  at  'JTRC  by  Briley  and  McDonald  (Ref.  9)  for  solution 
of  the  three-dimensional,  compressible  Navier-Stokes  equations.  This  procedure 
since  has  been  used  for  further  studies  of  rectangular  duct  flow  by  Briley, 
McDonald  and  Gibeling  (Ref.  10)  and  for  three-dimensional  combustor  flow  cal- 
culations by  Gibeling,  McDonald  and  Briley  (Ref.  18).  The  Navier-Stokes  por- 
tion of  the  present  study  is  based  upon  a two-dimensional  version  of  the  MINT 
code.  In  brief,  this  procedure  first  linearizes  the  equations  by  expanding  the 
solution  at  a known  time  level,  n,  to  represent  the  solution  at  time  level  n + 1. 
The  resulting  linear  equations  are  solved  using  an  alternating  direction 
implicit  (ADI)  method.  A description  of  this  numerical  procedure  is  presented 
in  Appendix  D.  More  detailed  descriptions  are  presented  in  Refs.  9i  1°  and  18. 


Results  of  Navier-Stokes  Computations 

Several  sample  cases  were  considered  in  assessing  the  ability  of  the  MINT 
code  to  make  predictions  of  the  shock  wave -boundary  layer  interaction  flow 
field.  The  shock  impingement  problem  was  considered  herein  since  some 
experimental  data  are  available  (Law,  Ref.  26),  and  this  problem  has  been 
considered  previously  by  other  investigators  (e.g.,  Refs.  29  and  48).  The  flow 
conditions  upstream  of  the  incident  shock  wave  were  chosen  to  correspond  as 
closely  as  possible  to  Law's  experiment  (Ref.  26).  Incident  shock  waves  of 
strength  p2/pp  = 2 (case  1,  9S  = 27.4°)  and  p2/pp  = 3 (case  2,  0C  = 33.8°) 
were  considered,  where  p^  is  the  static  pressure  upstream  of  the  shock,  p2  is 
the  static  pressure  downstream  of  the  shock,  and  9 s is  the  shock  angle.  The 
former  case  corresponds  approximately  to  the  shock  generator  angle  of  9.87° 
examined  by  Refs . 26  and  48 . 
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The  flow  field  in  the  present  analysis  is  divided  into  two  overlapping 
computational  domains  following  the  procedure  used  in  Ref.  48.  In  the  first 
domain  the  turbulent  boundary  layer  development  over  a flat  plate  was  calcu- 
lated. A 36x40  computational  grid  covering  a physical  domain  of  approximately 
56x606  was  employed  in  the  boundary  layer  calculation.  The  interaction  region 
is  contained  in  the  second  computational  domain.  For  case  1 a 36x61  computa- 
tional grid  was  employed  for  a domain  of  about  56x206,  and  for  case  2 a 
36x56  grid  was  employed  for  a domain  of  about  56x176 • A Roberts  grid  trans- 
formation was  used  in  the  normal  coordinate  direction  with  23  points  in  the 
inner  region  (Aym^n  = O.OO63  mm),  and  13  points  in  the  constant  mesh-spacing 
outer  region  (Ay  = 1.09  mm).  The  streamwise  grid  spacing  in  the  boundary 
layer  calculation  (first  computational  domain)  was  Ax  = 5.6  mm,  while  in  the 
interaction  domain  the  spacing  was  Ax  = 1.22  mm  for  case  1 and  Ax  = l.l4  mm 
for  case  2. 

In  this  feasibility  study,  computation  times  on  a CDC  6600  were  approxi- 
mately one  hour  per  case  to  reach  convergence.  However,  the  cases  were  run 
with  a generalized  code  capable  of  performing  calculations  both  in  three 
dimensions  and  with  orthogonal  curvilinear  coordinates  and  no  serious  effort 
was  made  to  optimize  either  the  code  or  the  time  step  for  this  particular 
problem.  If  the  code  were  rewritten  to  apply  specifically  to  a Cartesian 
system,  it  is  estimated  that  the  run  times  would  be  reduced  by  approximately 
a factor  of  two.  In  regard  to  optimization  of  the  time  step,  in  the  limited 
amount  of  time  and  effort  available  a detailed  study  concentrating  on  time 
step  optimization  for  this  type  of  problem  was  not  possible  and  a rather 
conservative  time  step  selected  which  doubtless  did  not  take  full  advantage 
of  the  large  time  step  capabilities  inherent  in  the  present  implicit 
procedure.  Finally,  solution  oscillations  related  both  to  the  shock 
representation  and  the  turbulence  model  (which  are  discussed  subsequently) 
may  also  have  increased  the  case  running  time.  Based  upon  these  initial  runs, 
it  is  estimated  that  a streamlined  Cartesian  deck  with  an  optimum  time  step 
could  perform  an  interaction  calculation  in  15  to  25  minutes  of  CDC  6600  run 
time  with  about  2200  mesh  points.  Finally,  it  should  be  noted  that  the  present 
calculations  were  run  with  a grid  definition  in  the  near  wall  region  much 
better  than  that  usually  used  in  stability  restricted  explicit  calculations. 

For  example,  the  grid  spacing  in  the  immediate  vicinity  of  the  wall  used  in 
the  present  effort  was  approximately  one-fifth  of  that  used  by  Shang,  Hankey, 
and  Law  (Ref.  48).  This  improvement  in  grid  definition  obtained  in  the 
present  effort  without  an  accompanying  penalty  in  computer  run  time  was  made 
possible  by  the  favorable  stability  properties  of  implicit  methods.  In  the 
upstream  boundary  layer,  the  grid  employed  in  the  present  effort  contained 
two  points  in  the  region  y+<  7.  Based  upon  our  experience  with  classical 
boundary  layer  solutions,  one  point  is  required  int.he  region  y+ < 7 to  obtain 
skin  friction  predictions  which  are  accurate  to  within  10  percent. 


The  sample  cases  pointed  out  several  deficiencies  in  the  "shock  capturing" 
treatment  of  shock  waves.  As  noted  by  Shang,  Hankey,  and  Law  (Ref.  44)  there 
is  significant  smearing  of  the  incident  shock  wave  in  the  relatively  coarse 
mesh  which  must  be  employed  in  a Navier-Stokes  calculation.  This  shock 
resolution  problem  appears  to  be  more  pronounced  in  the  current  implicit 
procedure  than  it  was  in  explicit  procedures  based  upon  MacCormack's  method. 
However,  this  is  not  surprising  since  MacCormack's  method  was  developed  for 
shock  calculations  over  a period  of  many  years,  whereas  the  present  work  is 
the  first  application  of  the  MINT  procedure  to  flow  fields  containing  shock 
waves.  The  present  results  should  not  be  construed  as  indicating  an  inherent 
limitation  of  implicit  schemes  to  resolve  shock  waves,  they  imply  that  at 
th'  • stage  of  development  they  do  not  capture  shocks  quite  as  well  as  some 
of  the  explicit  schemes.  Obviously,  much  further  work  aimed  at  applying 
implicit  techniques  to  shock  wave  problems  needs  to  be  done . 

Plots  of  skin  friction  coefficient  and  surface  pressure  are  shown  in 
Figs.  20a  and  20o  for  case  1 (pp/Pq  = 2)  and  in  Figs.  21a  and  21b  for  case  2 
(P2/P1  = 3)*  Although  both  experimental  data  as  well  as  the  explicit  calcu- 
lations of  Shang,  Hankey,  and  Law  (Ref.  48)  indicate  separation  under  these 
conditions,  the  present  calculations  do  not.  There  are  several  possible 
reasons  for  this  discrepancy.  First,  as  previously  stated,  the  shock  smearing 
problem  in  the  present  implicit  calculation  appears  to  be  more  severe  than  in 
explicit  calculations  performed  to  date.  Secondly,  a more  sophisticated 
turbulence  model  than  that  which  was  used  may  be  required.  Finally,  in 
comparison  with  the  calculations  of  Ref.  48,  the  present  procedure  had  a much 
better  definition  of  the  flow  in  the  wall  region.  It  is  not  clear  what  the 
effect  of  improving  the  grid  definition  in  the  explicit  method  would  have  upon 
the  explicit  calculation  results  and,  therefore,  a direct  comparison  between 
the  present  implicit  calculation  and  the  explicit  calculation  of  Ref.  48  is 
somewhat  difficult.  The  experimental  surface  pressure  distribution  for  case  1 
is  also  shown  in  Fig.  20b,  from  which  it  is  apparent  that  the  initial  rapid 
rise  of  surface  pressure  which  accompanies  separation  is  not  predicted  very 
well  by  the  present  calculations.  This  seems  to  be  directly  related  to  the 
smearing  of  the  incident  shock  wave  and  is  beiieved  to  be  closely  connected 
with  the  lack  of  separation  in  this  case. 

Since  the  pressure  ratio  p2/p^  = 2 case  did  not  separate,  a stronger 
shock  calculation  having  a pressure  ratio  of  3 was  considered.  As  shown  in 
Fig.  21b,  the  shape  of  the  predicted  pressure  rise  showed  the  rapid  initial 
increase  expected  in  interaction  flow  fields.  Furthermore,  the  level  of  the 
downstream  wall  pressure  was  essentially  the  inviscid  value.  The  predicted 
length  of  the  separated  region  was  shorter  than  expected,  a result 
consistent  with  the  lack  of  separation  observed  in  the  pressure  ratio  po/Pi  = 2 


calculation.  Finally,  the  skin  friction  distribution  in  the  separated  case 
(Fig.  21a,  P2/P1  = 3)  is  also  qualitatively  reasonable,  although  the  increase 
in  skin  friction  downstream  of  the  incident  shock  was  unexpected.  This  may  be 
due  to  the  turbulence  model  incorporated  in  the  present  procedure. 

Several  problems  have  been  encountered  in  the  present  study  which  require 
further  investigation.  First  of  all,  shock  resolution  in  the  implicit  frame- 
work must  be  improved,  since  the  resolution  obtained  with  the  present  formula- 
tion is  apparently  not  as  good  as  that  obtained  using  MacCormack's  alternating  - 
direction-explicit  method  (e.g.,  Refs.  29  and  48).  Also,  some  oscillation  of 
the  solution  about  the  steady-state  has  been  observed,  especially  in  the  sub- 
layer. This  problem  may  be  partially  related  to  both  the  turbulence  model  and 
the  choice  of  time  step.  Further  investigation  into  the  optimal  choice  of 
the  time  step  for  this  implicit  procedure  should  be  carried  out,  since  it  is 
believed  that  a significant  improvement  in  the  convergence  rate  can  be 
realized.  Finally,  the  present  calculations  were  obtained  using  primitive  o 
variables  ( D,  u,  v,  h)  rather  than  conserved  variables  (p,  pu,  pv,  p(e  + fq  ) ) 
which  may  have  resulted  in  additional  smearing  of  the  incident  shock  wave. 
Hence,  in  future  studies  consideration  should  be  given  to  using  conserved 
variables . 
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Figure  20.  Navier  — Stokes  Interaction  Calculation,  P2/P1  = 2 
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SECTION  III 
CONCLUSIONS 


BOUNDARY  LAYER  ANALYSIS 

'The  strong  interaction  boundary  layer  approach  has  been  applied  to  the 
shock  wave  turbulent  boundary  layer  interaction  problem.  A well  proven 
weak  interaction  boundary  layer  procedure  was  modified  to  allow  for  strong 
interaction  solutions.  The  modifications  included  the  incorporation  of  an 
outer  flow  analysis  and  the  revision  of  the  governing  equations  to  allow  for 
flow  in  separated  regions.  When  the  turbulence  model  originally  developed 
for  weak  interaction  flows  was  used,  poor  agreement  was  found  between  the 
calculations  and  data.  This  is  consistent  with  the  results  reported  by  Shang 
and  Hankey  (Ref . 47).  Based  upon  computational  turbulence  studies  and  recent 
experimental  data,  a modification  to  the  basic  turbulence  model  was  developed 
for  boundary  layers  containing  recirculating  flow.  Using  this  modified  tur- 
bulence model,  good  agreement  was  found  between  the  calculations  and  test 
data. 


NAVIER -STOKES  ANALYSIS 

An  implicit  finite-difference  Navier-Stokes  analysis  has  been  applied  to 
the  shock  wave -boundary  layer  interaction  problem.  Despite  the  fact  that  the 
code  was  written  in  general  curvilinear  coordinates  and  was  not  streamlined 
for  the  interaction  problem,  solutions  were  obtained  in  reasonable  run  times. 
Indications  are  that  streamlining  the  code  and  optimizing  the  time  step  could 
decrease  the  run  time  by  a factor  of  four.  Although  solutions  for  strong 
incident  shocks  showed  a qualitatively  correct  behavior  (i.e.,  skin  friction 
and  wall  pressure  distributions  were  in  qualitative  agreement  with  data),  the 
procedure  severely  smoothed  the  incident  shock  wave  thus  suppressing 
separation  for  moderate  shocks  and  underpredicting  the  extent  of  separation 
for  strong  snocks.  Therefore,  calculation  of  shock  waves  in  an  implicit 
method  requires  further  investigation.  In  addition,  further  work  in  regard  to 
turbulence  modeling  and  optimum  time  step  specification  is  warranted. 
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APPENDIX  A 

WEAK  INTERACTION  CALCULATION  PROCEDURE 

The  procedure  used  to  solve  the  momentum  and  energy  equations  is  a 
Hartree-Womersley  calculation  procedure  similar  to  that  described  by  Smith 
and  Clutter  (Ref.  50) . In  the  calculation  procedure  the  differential  equa- 
tions are  transformed  to  a form  more  convenient  for  computer  solution  by 
introducing  new  variables 
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where  the  primes  indicate  differentiation  with  respect  to  T]  , T°  is  the 
stagnation  temperature,  Tp>gp  is  a specified  constant  reference  temperature, 

and  3+  is  a length  scale  usually  taken  to  be  the  displacement  thickness.  In 
certain  flows,  such  as  the  boundary  layer  developing  on  the  cold  wall  nozzle, 
the  displacement  thickness  may  become  small  or  even  negative  and  in  these  cases 
6+  is  taken  to  be  a linear  combination  of  the  displacement  thickness  and  a 
constant  reference  length  such  that  6f  remains  positive. 

The  momentum  equation,  Eq.  (9) , and  the  energy  equation,  Eq.  (10),  are 
solved  by  first  eliminating  pv  through  the  continuity  equation,  Eq.  (3), 
and  replacing  the  variables  T°,  p,  and  u by  G , F,  and  fl.  The  streamwise 
derivatives  are  then  replaced  by  finite  differences  leading  to  equations  of 
the  form 

A3  F'"+  A2  F"+  A,  F'+  A0  F = A4  (A-5) 

B3  G'"  + BjG  + B/  G*  : B4  (A-o) 

I 

where  An  and  Bn  are  functions  of  F,  G and  their  derivatives.  The  equations 
are  linearized  by  assuming  values  for  F,  G1  and  their  derivatives  based  upon 
the  solution  at  the  previous  streamwise  stations  and  the  resulting  linear  dif- 
ferential equations  are  solved  by  Gaussian  elimination  using  the  boundary 
conditions 
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F '(S)  = 0 
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where  6th  the  thickness  of  the  thermal  boundary  layer.  In  the  calculations 
presented  in  the  present  report,  uw  was  always  taken  as  zero.  Having  obtained 
the  solution  of  the  linearized  equations,  the  distributions  of  F and  G’  are  com- 
pared to  the  distributions  used  to  evaluate  the  nonlinear  coefficients,  An  and 
Bn.  If  the  old  and  new  distributions  agree  to  within  a specified  tolerance, 
the  procedure  moves  to  the  next  streamwise  station.  If  the  two  do  not  agree, 
the  procedure  is  repeated.  In  the  case  of  turbulent  flow,  the  coefficients 
An  and  Bn,  also  depend  on  the  turbulent  kinematic  viscosity,  vp,  and  a turbu- 
elnce  kinetic  energy  equation  is  included  in  the  iteration  loop. 


APPENDIX  B 


BOUNDARY  LAYER  MODIFICATIONS  FOR  6*  SPECIFIED 

To  specify  the  displacement  thickness,  o*s  and  evaluate  the  static 
pressure  from  the  momentum  equation,  a relation  containing  6*  was  derived  to 
replace  the  pressure  gradient  terms  in  the  momentum  equation.  A parameter, 
Q,  was  introduced  to  represent  the  naturally  occurring  term: 


Q = 


(B-l ) 


Q was  then  treated  as  one  of  the  dependent  variables  and  terms  containing 
products  of  Q and  other  dependent  variables  were  linearized  in  the  usual 
manner. 


From  isentropic  relations  it  can  be  shown  that 
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Expanding  the  definition  of  Q yields 
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which  is  of  the  linear  form: 
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When  specifying  the  displacement  thickness,  §*>  Tj^  and  T2  in  Eq.  (B-5)  are 
evaluated  and  the  linear  relation  in  Q is  used  to  replace  Bue/dx  in  the  pres- 
sure gradient  terms  of  the  momentum  equation. 
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The  momentum  equation  is  solved  for  the  6*  specified  case  by  linearizing 
the  partial  differential  equations  and  replacing  derivatives  by  differences 
at  each  grid  point.  The  resulting  set  of  algebraic  equations  is  then  solved 
for  the  value  of  the  stream  function,  F,  at  each  grid  point  and  for  the 
parameter  Q.  The  solution  is  obtained  by  a modified  Gaussian  elimination 
which  solves  a quin-diagonal  matrix  with  an  additional  column. 


APPENDIX  C 


ORIGINAL  TURBULENCE  MODEL  USING  THE  TURBULENCE  KINETIC  ENERGY  EQUATION 


As  shown  in  Ref.  46,  the  boundary  layer  approximation  to  the  turbulence 
kinetic  energy  equation  is  given  by 
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All  calculations  reported  in  the  investigation  were  made  with  the  usual 
assumption  of  zero  pressure-dilitation  contribution  to  the  energy  balance 
(Ref.  8)  unless  otherwise  stated.  The  turbulence  model  is  developed  by 
integrating  Eq.  (C-l)with  respect  to  y between  the  limits  y = 0 and  y = 6 
which  leads  to 
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where 


E = [ -j-  q2  (p  u - p v)  - p'v'  (pv) ' q 2]e  (c 

Following  Townsend  (Ref.  48)  and  Bradshaw  and  Ferris  (Ref.  8)  structural 
coefficients  an  and  L are  introduced,  together  with  a mixing  length  £; 
these  scales  are  defined  as 

- uV  = a,  q2  , u'2:a2q7,  v'2  ^ a3  q 
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For  fully-developed  turbulence  the  structural  coefficients  a-^,  a0,  and  a^ 
are  assumed  constant  having  values  0.15,  0.50,  and  0.20,  respectively  (Refs. 
7 and  31).  Using  Eq.  (C-4),Eq.  (C-3)  is  Put  in  the  form 


where 
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where  H is  a nondimens ional  transverse  distance  y/6+,  6+  is  an  arbitrary 
reference  length,  and  6 the  boundary  layer  thickness. 

The  left-hand  side  of  Eq.  (C-5)  represents  the  streamwise  rate  of  change 
of  turbulence  kinetic  energy  and  is  derived  directly  from  the  turbulence 
kinetic  energy  advection  term.  The  term  Peu^2  represents  the  integral  of 
turbulence  production  minus  dissipation  and  Peu3^^  is  the  normal  stress  pro- 
duction. The  terms  designated  by  E are  turbulent  source  terms  resulting  from 
disturbances  imposed  upon  the  boundary  layer  by  the  free  stream.  As  shown  in 
Eq.  (C-3),E  is  the  sum  of  two  major  contributions,  the  first  (q^/2)  ( Pu^/$x-pv) 
representing  the  free-stream  velocity  disturbance  ^i.e.,  free-stream  turbulence 
entrained  by  the  boundary  layer)  and  the  second,  P v'  + (pv)'q^/2,  representing 
the  direct  absorption  of  acoustic  energy. 

For  fully-developed  turbulent  flow,  as  in  Ref.  31,  L and  £ are  given 
by 

-g-  = 0.1  tanh  [<cy  /( 0 I 8)]  (c-9) 

-g-  = tonh  [*y/ 4»]  (c-io) 

where  is  the  "wake"  value  of  the  mixing  length  at  any  particular  streamwise 
station.  Although  Eqs . (C-9)and  (C-10)  give  accurate  representations  of  £ and 
L through  most  of  the  turbulent  boundary  layer,  it  is  well-known  that  they 
overestimate  the  length  scales  within  the  viscous  sublayer  and  are  somewhat 
inaccurate  at  low  Reynolds  numbers.  Following  McDonald  and  Fish  (Ref.  33) 
the  experimentally  observed  damping  effect  in  the  viscous  sublayer  is  modeled 
by  assuming  intermittent  turbulence  within  the  sublayer  leading  to  the  relation 
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In.  Eq.  (C-6)  r is  the  intermittency  factor,  the  damping  factor,  and  the 
subscript  T indicates  the  value  with  turbulent  flow.  Obviously,  lb  is 
equal  to  the  square  root  of  p.  As  in  Ref.  33,  the  present  investigation 
assumes  that  the  damping  distributes  normally  about  a mean  height  y+(y+= 
with  a standard  deviation  leading  to  the  equation 

- p 1/2  {(y*  -r  )/*}  <c-12) 

where  P is  the  normal  probability  function;  y^  is  taken  as  23,  and  o as  8. 

A detailed  discussion  of  the  sublayer  damping  treatment  is  presented  in  Ref. 

33.  In  the  present  calculations  the  von  Karman  constant  k was  taken  to  be 

0.43. 


In  regard  to  the  low  Reynolds  number  effects,  Coles  (Ref.  12)  has  observed 
and  correlated  the  departure  of  the  mean  velocity  profile  of  a flat  plate 
turbulent  boundary  layer  from  the  usual  similarity  laws  known  to  hold  at 
higher  Reynolds  numbers.  Using  Coles'  correlation  of  the  mean  velocity  pro- 
file in  the  low  Reynolds  number  regime,  McDonald  (Ref.  31)  integrated  the 
boundary  layer  equations  of  mean  motion  to  obtain  local  distributions  of 
turbulent  shear  stress  and  evaluated  the  local  mixing  length  distributions 
from  the  assumed  mean  velocity  distribution  and  the  computed  shear  stress 
distributions.  Based  upon  these  calculations,  a low  Reynolds  number  correc- 
tion for  the  dissipation  length  of  the  form 

L = !_<*,[  I + e*P(-l  63  In  Rg+  9 7)]  (C-13) 

was  derived  where  Loo/6  is  given  by  Eq.  (C-9).  In  the  calculations  presented 
in  the  present  report  the  dissipation  length  used  was  obtained  by  multiplying 
Eq.  (C -13)  by  the  sublayer  damping  factor , lb. 

When  numerical  values  of  the  structural  coefficients  an  are  specified, 

Eqs , (c-10) . (C-12 ) , and  (C-13)  are  used  to  represent  L and  i , and  the  pressure 
dilitation  is  either  neglected  or  modeled,  the  turbulence  kinetic  energy 
equation,  Eq.  (C-f ),  becomes  an  ordinary  differential  equation  with  the  depen- 
dent parameter  &°(x)  which  is  solved  in  conjunction  with  the  boundary  layer 
momentum  and  energy  equations  to  predict  the  development  of  both  the  mean  flow 
field  and  the  turbulent  shear  stress. 

In  addition  to  including  the  turbulence  kinetic  energy  equation  in  the 
set  of  equations  governing  the  boundary  layer  development  it  is  necessary 
to  specify  a model  for  turbulent  heat  flux  contribution,  "pCpv'T'.  As  pre- 
viously stated,  I si  the  present  procedure,  v'T'  is  specified  by  assuming  a 
turbulent  l’randtl  number,  Fr^,  which  relates  the  velocity- temperature  corre- 
lation, v'T',  to  the  Reynolds  stress,  u'v',  through  Eq.  (8).  The  turbulent 
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Prandtl  number  distribution  used  in  the  present  procedure  varies  with  distance 
from  the  wall  in  the  manner  suggested  by  Meier  and  Rotta  (Ref.  35 ) • As  this 
juncture  it  should  be  pointed  out  that  an  alternative  procedure  can  be  used 
to  determine  v'T',  based  upon  an  easily  derived  conservation  equation  for 
either  the  quantity  T'^  or  the  correlation,  v'T' , which  is  similar  in  form 
to  the  turbulence  kinetic  energy  equation,  Eq.  (23).  However,  to  solve 
this  new  conservation  equation  it  is  necessary  to  assume  a universal  struc- 
ture relating  quantities  analogous  to  dissipation,  production,  etc.  While 
sufficient  experimental  data  exists  to  allow  valid  modeling  of  the  required 
terms  for  the  turbulence  kinetic  energy  equation,  the  existing  data  does  not 
indicate  how  proper  modeling  could  be  carried  out  for  the  v'T'  conservation 
equation.  Thus,  at  least  for  the  present,  the  approach  based  upon  a turbu- 
lent Prandtl  number  appears  preferable  to  an  approach  based  upon  the  v'T'  con- 
servation equation. 


APPENDIX  D 


THE  MINT  CALCULATION  PROCEDURE 


The  MINT  calculation  procedure  represents  an  efficient  solution  to  the 
multi-dimensional,  time -dependent  Navier-Stokes  equations.  The  procedure  is 
presented  in  detail  in  Refs.  9»  10  lb;  for  convenience  it  is  presented  in 
condensed  form  in  this  appendix.  As  a prelude  to  this  discussion  it  is  con- 
venient to  define  a difference  notation.  For  the  present  it  is  assumed  that 
the  flow  region  is  three-dimensional  and  is  discretized  by  grid  points  having 
equal  spacings  A x^,  A X2  and  A x^  in  the  x^,  x2  and  x^  directions  respectively; 
in  addition  the  time  step  is  At.  Reduction  to  two  dimensions  will  be  obvious 
as  the  discussion  proceeds  and  provisions  for  nonuniform  grids  will  be  intro- 
duced at  the  end  of  this  appendix. 

The  subscripts  i,  j,  k and  superscript  n are  grid  point  indices  associated 
with  x-,,  x2,  and  x^,  and  t,  respectively.  Thus  9^  j k denotes  0 (xii>  x2j> 
xg^,  t*)  where  0 can  represent  any  of  the  dependent  variables . The  subscripts 
are  frequently  omitted  if  clarity  is  preserved,  so  that  $ is  equivalent  to 
^i  j k.  For  convenience,  the  following  shorthand  difference-operator  notation 
is  used  for  derivative  difference  formulas : 


2 A x 


dx 


0(Ax| 


(D-l) 
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with  analogous  definitions  for  62,  60,  and  5^.  It  is  assumed  that  the 

solution  is  known  at  the  n level,  tn,  and  is  desired  at  the  (n+l)  level, 

, n+1 


LINEARIZATION 


The  large  time-step  capabilities  of  implicit  methods  can  place  great 
demands  on  the  linearization  technique  employed.  Indeed,  the  favorable 
stability  properties  of  implicit  methods  can  be  severly  compromised  by  an 
inadequate  linearization.  The  technique  used  in  the  MINT  code  permits  the 
implicit  solution  of  coupled  nonlinear  equations  in  one  space  dimension  by 
a one-step  noniterative  procedure.  This  feature  is  retained  for  multidimen- 
sional problems  by  using  ADI  techniques.  The  linearization  is  accurate  when 
variables  change  by  relatively  small  amounts  during  a time  step,  and  conse- 
quently, the  accuracy  of  the  linearization  can  be  controlled  by  varying  the 
time  step.  The  linearization  technique  is  also  convenient  for  the  implicit 
treatment  of  coupled  nonlinear  boundary  conditions,  and  this  latter  feature 
has  been  found  to  have  a highly  favorable  effect  on  the  stability  of  the  over- 
all method  (Ref.  9)« 

For  demonstration  purposes  the  technique  now  will  be  described  for  the 
fallowing  first  order  equation  in  one  spatial  variable,  0 (x,t) 


d <p 

JT 


= F (<£ ) J-  G (<*>) 


(D-3) 


The  procedure  is  based  on  an  expansion  of  nonlinear  implicit  terms  about 
the  known  time  level,  tn,  and  leads  to  a one-step,  two-level  scheme  which, 
being  linear  in  unknown  (implicit ) quantities,  can  be  solved  efficiently  without 
iteration.  The  technique  is  easily  extended  to  treat  coupled  systems  of 
equations  and  second-order  spatial  derivatives.  The  difference  approximation 
is  derived  from  the  following  backward  time-difference  replacement  of  Eq. 

(D-3) 


n ♦ i 

F (4>)  J- G (<t>)  + 0 (At)  (D-U) 

dx 

i 

where  the  spatial  differencing  of  the  bracketed  term  is  as  yet  unspecified. 
Making  use  of  chain-rule  differentiation,  the  bracketed  term  in  Eq.  (D-4 ) is 
expanded  about  t ; the  result  is  then  differenced  using  forward  time  differ- 
ences and  centered  spatial  differences  to  obtain  the  following  implicit  differ- 
ence scheme: 
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On  examination,  it  can  be  seen  that  Eq.  (D-5)  is  linear  in  and  that 

all  other  quantities  are  either  known  or  evaluated  at  the  n level.  Because 
of  the  spatial  difference  operator,  6X,  Eq.  (D-5)  contains  ^+1>  and 

consequently,  the  system  of  linear  equations  generated  by  writing  Eq. 
(D-5)  at  each  of  the  grid  points,  x^,  must  be  solved  simultaneously  as  an 
implicit  system.  The  implicit  system  of  equations  can  be  written  in 
tridiagonal  matrix  form,  and  can  therefore  can  be  solved  easily  and  efficiently 
by  standard  techniques  for  tridiagonal  systems  (see,  e.g.,  Ref.  23).  The 
tridiagonal  matrix  structure  emerges  from  writing  Eq.  (D-5)  in  the  following 
form: 


n , n»l  , n ,n*i  n n*i  n 

°i  *l-i  +b,  *i  +Ci  *,+,=di 


(D-6) 


where  the  coefficients  contain  only  n-level  quantities.  When  applied  at 
successive  grid  points,  Eq.  (D-6)  generates  a tridiagonal  system  of  equations 
for  f+1. 


APPLICATION  OF  THE  METHOD 


The  extension  of  the  numerical  method  to  more  than  one  spatial  dimension 
is  based  upon  an  alternating  direction  implicit  (ADI)  technique.  The  technique 
is  an  application  of  the  general  procedure  developed  by  Douglas  and  Gunn 
(Ref.  15)  in  which  the  linearization  technique  described  previously  is  applied 
to  the  coupled  system  of  governing  equations,  Eqs.  (2b)  - (26). 


These  equations  are  written  in  backward  time  difference  form,  and  non- 
linear implicit  terms  are  linearized  by  expansion  about  tn0.  The  viscous 
force  terms  in  Eq.  (25)  which  contain  mixed  derivatives  (i.e.,  A /Ax  . Ax ^ for 
i / j)  are  most  easily  treat  explicitly  by  evaluation  at  time  level  n. 
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Although  mixed  derivatives  can  be  differenced  implicitly  within  the  Douglas-Gunn 
framework,  this  would  increase  the  number  of  intermediate  steps  and  thereby 
complicate  the  solution  procedure.  Previous  experience  with  the  method  in 
Refs.  9,  10  and  18  indicates  that  the  explicit  treatment  of  these  viscous 
terms  had  no  observable  effect  on  stability. 

The  difference  equations  obtained  by  the  procedure  outlined  above  repre- 
sent a linearized  backward  difference  scheme.  The  equations  can  be  arranged 
according  to  time  and  space  derivatives,  and  written  in  the  following  matrix 
operator  form  (Ref.  9): 


=n  ( 4>  - \ = n - n »i  = n Tn*'  » n tD*i  _n 

A D,  <p  + D2  4>  + D3  4>  + s" 


(D-7) 


=n 

Here  A is  a (mxm)  matrix  containing  the  time  derivative  coefficients,  where 
m is  the  number  of  equations_j^eing  solved;  0 is  the  column  vector  of  the  depen- 
dent variables;  Dj_,  D2,  and  D3  are  (mxm)  matrices  containing  three-point  dif- 
ference operators  associated  with  the  coordinate  directions  x-^,  x,J5  and  x.,, 
respectively;  and  S is  a column  vector  containing  only  n-level  terms.  Since 
the  multidimensional  implicit  system  with  coefficients  generated  by  Eq.  (D-7) 
is  difficult  to  solve,  the  Douglas-Gunn  (Ref.  15)  technique  is  applied  to 
Eq.  (D-7)  to  generate  an  ADI  scheme.  With  the  observation  that  the  Douglas- 
Gunn  procedure  is  being  applied  to  a coupled  system  of  equations,  the  follow- 
ing three-step  scheme  is  obtained.  (For  two  spatial  dimensions  the  technique 
collapses  into  a two-step  scheme.) 


= n 7 * = n , n =.  nr-n  — n 

= D,  9 + D2  9 + D3  9 + S 


(D-8) 


A 


(D-9) 


-n 

A 


= n 1' 
D3^> 


+ s 


(D-10) 


a 


JHt’  ■.-JHt'Jt 

where  0,0,  and  0 are  the  intermediate  solutions.  Note  that  at  each 
step  of  the  scheme,  one  more  coordinate  direction  is  treated  implicitly,  and 
that  the  most  recent  approximation  to  0 is  not  always  used,  as  this  would 
adversely  affect  the  stability. 

The  effort  involved  in  the  actual  programming  and  solution  of  Eqs.  (D-8)  - 
(D-10)  is  greatly  reduced,  and  the  computer  storage  requirements  are  halved 
by  subtracting  Eq.  (D-8)  from  Eq.  (D-9)  and  Eq.  (D-9)  from  Eq.  (D-10),  so 
that  Eqs.  (D-9)  and  (D-10)  have  the  following  simplified  form: 

A"(il£l)=B"2r-.n  own 


For^three  spatial  dimensions  0 represents  the  solution  at  time  (n+1), 

0 ; for  two  spatial  dimensions  Eq.  (D-12)  is  eliminated  and  0 * represents 

the  solution. 

On  examination,  it  can  be  seen  that  the  difference  equations  (D-8),  (D-ll), 
(D-12 ) are  linear  in  the  * -level  quantities.  At  the  kth  step  in  the  procedure 
there  are  m equations  at  each  gjid  point  (x-^,  x2,  x^);  because  of  the  spatial 
difference  operators  (6^  and  6^  ) these  equations  contain  the  dependent  vari- 
ables at  x^  and  at  each  of  the  two  adjacent  grid  points  in  the  x^-direction. 
Consequently,  the  difference  equations  nmst  be  solved  as  an  implicit  system. 

It  should  be  recognized  that  upon  application  at  a successive  number  of  grid 
.uation  general^  a blog^.-tridi  agonal  system  of  algebraic 
(i.e.,  0,0  , or  0 FOR  K = 1,  2,  3).  After  appro- 

priate treatment  of  boundary  conditions,  each  system  can  be  solved  efficiently 
using  a standard  block  elimination  method  such  as  the  matrix  factorization 
method.  The  method  used  in  the  present  study  is  closely  related  to  the 
matrix  factorization  method  and  simply  consists  of  Gaussian  elimination  for 
a tridiagonal  matrix,  but  whore  the  elements  of  the  tridiagonal  matrix  are 
(mxm)  submatricos  rather  than  scalars,  where  m is  the  number  of  governing 
equations . 


points,  x^,  each, 
equations  for  0 ' 


-'7 
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The  solution  procedure  for  a single  time  step  is  as  follows; 

(1)  During  the  first  step  of  the  ADI  procedure,  Eq.  (D-8)  is  applied  at 

successive  x^-direction  rows  of  grid  points  to  provide  one-dimensional 
implicit  systems_of  equations.  These  systems  are  generated  by  the 
operator  (A  /At-Bp).  The  implicit  systems  can  be  arranged  in  block- 
tridiagonal  form  and  solved  as  indicated  previously.  Since  there 
are  m governing  equations,  the  block-tridiagonal  systems  have 
(mxm)  square  matrices  as  the  block  elements. 


Although  this  discussion  has  assumed  the  grids  to  be  uniform,  extension 
to  nonuniform  grids  is  straightforward.  When  a Roberts  transformation  is 
used  (see  Eqs.  (4l)  - (4'd)),  the  first  and  second  derivatives  for  unequally 
spaced  grid  points  still  le -d  to  three  point  difference  operators.  Thus  the 
resulting  form  of  the  matrices  is  unchanged  and  the  solution  procedure  is 
identical  to  that  used  for  the  equally  spaced  case. 
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